Simulating electron motion - differential equation with adaptive step size in python

Question

I am trying to simulate how the oscillating electric field of an intense laser will push around an electron that is near the Coulomb potential of a +1 ion. The laser field is

E = Eo sin(wt), in the y direction.

and the Coulomb potential is

F = ke q1*q2/r^2, in the r direction.

The strong electric field causes the electron to tunnel ionize, so the initial condition of the electron is to be displaced from the atom in the y-direction. Then, the electron is pushed back and forth by the laser field and has a chance to interact with the Coulomb potential. I want to simulate how the Coulomb potential affects the flight of the electron. The simulations need to be in three dimensions, because I eventually want to include more complex laser fields that push the electron in the x and y directions and the electron can start with momentum in the z direction.

At first, I thought that this would be easy. Below is the code that I used to step through time in very small steps (1e-18 sec). When the electron is not near the ion, this works fine. However, for electrons that pass close to the ion, the results depend strongly on the time-step used in the simulations. If I make the time-step smaller, the calculations take a very long time.

So, I think in this case I am supposed to use an adaptive timestep. Also, from what I have read, the Runge-Kutta methods are supposed to be superior to the simple approach I am using. However, I don't think that scipy.odeint applies to three-dimensions. Any ideas on how to improve the accuracy and speed of these simulations?

Here is the figure showing how the time-step has a huge impact on the results (a bad thing):

And here is my code:

import numpy as np
import matplotlib.pyplot as plt

q = 1.602e-19       #Coulombs   Charge of electron
h_bar = 1.054e-34   #J*s        Plank's Constant div by 2Pi
c = 3.0e8           #m/s        Speed of light
eo = 8.8541e-12     #C^2/(Nm^2) Permittivity of vacuum
me = 9.109e-31      #kg         Mass of electron
ke = 8.985551e9     #N m^2 C-2  Coulomb's constant

def fly_trajectory(wavelength,intensity,tb=0,pulseFWHM=40.e-15,
                   final_time=100e-15,time_step=.001e-15,Ip=15.13,v0=(2e4,0,0)):
    #Intensity is in w/cm2. Ip is in eV. Otherwise it's SI units throughout.
    #The electric field of the later is in the y-direction

    Ip = 15.13 * q  #Calculate the ionization potential of the atom in Joules
    Eo = np.sqrt(2*intensity*10**4/(c*8.85e-12)) # Electric field in V/m
    w = c/wavelength * 2. * np.pi #Angular frequency of the laser

    times = np.arange(tb,final_time,time_step)
    Ey = Eo*np.sin(w*times) * np.exp(-times**2/(2*(pulseFWHM / 2.35482)**2)) 
    Eb = Ey[0] #E-field at time of birth (time of tunneling)

    if Eb == 0: return 0,0 #No field --> no electrons
    tunnel_position = -Ip / (Eb*q)    
    x,y,z = 0,tunnel_position,0
    vx,vy,vz = v0

    y_list = np.zeros(len(times)) #store the y-values for later

    for index in range(0,len(times)):
        r = np.sqrt(x**2+y**2+z**2)
        rx = x/r; ry = y/r; rz=z/r

        Fcy = -q**2 * ke/(r**2) * ry
        ay = Ey[index]*(-q)/me + Fcy/me #only y includes the laser
        vy = vy + ay*time_step 
        y = y + vy * time_step

        Fcx = -q**2 * ke/(r**2) * rx
        ax = (-q)/me + Fcx/me
        vx = vx + ax*time_step  
        x = x + vx * time_step

        Fcz = -q**2 * ke/(r**2) * rz
        az = (-q)/me + Fcz/me
        vz = vz + az*time_step 
        z = z + vz * time_step

        y_list[index] = y

    return times,y_list

for tb in np.linspace(0.25*2.66e-15,0.5*2.66e-15,5):
    print tb
    times,ys = fly_trajectory(800e-9,2e14,tb=tb,time_step=.01e-15)
    plt.plot(times,ys,color='r')

    times,ys = fly_trajectory(800e-9,2e14,tb=tb,time_step=.001e-15)
    plt.plot(times,ys,color='b')


#plot legend and labels:
plt.plot(0,0,color='r',label='10e-18 sec step')
plt.plot(0,0,color='b',label='1e-18 sec step')
plt.xlim(0,10e-15); plt.ylim(-1e-8,1e-8)
leg = plt.legend(); leg.draw_frame(False)
plt.xlabel('Time (sec)')
plt.ylabel('Y-distance (meters)')
plt.show()

Answer 1

As Warren Weckesser suggested, I can simply follow the Scipy cookbook for the coupled mass-spring system. First, I need to write my "right side" equations as:

x'  = vx
y'  = vy
z'  = vz
vx' = Ac*x/r
vy' = Ac*y/r + q*E/m
vz' = Ac*z/r 

where Ac=keq^2/(mr^2) is the magnitude of the acceleration due to the Coulomb potential and E is the time-dependent electric field of the laser. Then, I can use scipy.integrate.odeint to find the solutions. This is faster and more reliable than the method that I was using previously.

Here is what the electron trajectories look like with odeint. Now none of them fly away crazily:

And here is the code:

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

q  = 1.602e-19    #Coulombs   Charge of electron
c  = 3.0e8        #m/s        Speed of light
eo = 8.8541e-12   #C^2/(Nm^2) Permittivity of vacuum
me = 9.109e-31    #kg         Mass of electron
ke = 8.985551e9   #N m^2 C-2  Coulomb's constant

def tunnel_position(tb,intensity,wavelength,pulseFWHM,Ip):
    Ip = 15.13 * q  
    Eb = E_laser(tb,intensity,wavelength,pulseFWHM) 
    return -Ip / (Eb*q) 

def E_laser(t,intensity,wavelength,pulseFWHM):
    w = c/wavelength * 2. * np.pi #Angular frequency of the laser
    Eo = np.sqrt(2*intensity*10**4/(c*8.85e-12)) # Electric field in V/m
    return Eo*np.sin(w*t) * np.exp(-t**2/(2*(pulseFWHM / 2.35482)**2))

def vectorfield(variables,t,params):
    x,y,z,vx,vy,vz = variables
    intensity,wavelength,pulseFWHM,tb = params
    r = np.sqrt(x**2+y**2+z**2)
    Ac = -ke*q**2/(r**2*me)
    return [vx,vy,vz,
         Ac*x/r,
         Ac*y/r + q/me * E_laser((t-tb),intensity,wavelength,pulseFWHM),
         Ac*z/r]

Ip  = 15.13   # Ionization potential of Argon eV
intensity  = 2e14
wavelength = 800e-9
pulseFWHM  = 40e-15

period  = wavelength/c
t = np.linspace(0,20*period,10000)

birth_times = np.linspace(0.01*period,0.999*period,50)
max_field   = np.max(np.abs(E_laser(birth_times,intensity,wavelength,pulseFWHM)))

for tb in birth_times:
    x0  = 0 
    y0  = tunnel_position(tb,intensity,wavelength,pulseFWHM,Ip)
    z0  = 0
    vx0 = 2e4
    vy0 = 0
    vz0 = 0

    p = [intensity,wavelength,pulseFWHM,tb]
    w0 = [x0,y0,z0,vx0,vy0,vz0]

    solution,info = scipy.integrate.odeint(vectorfield,w0,t, args=(p,),full_output=True)
    print 'Tb: %.2f fs - smallest step : %.05f attosec'%((tb*1e15),np.min(info['hu'])*1e18)

    y = solution[:,1]

    importance = (np.abs(E_laser(tb,intensity,wavelength,pulseFWHM))/max_field)
    plt.plot(t,y,alpha=importance*0.8,lw=1)

plt.xlabel('Time (sec)')
plt.ylabel('Y-distance (meters)')

plt.show()