Accurate numerical simulation of radiation reaction effects in strong electromagnetic fields

0.65 MB·

by  N. V. Elkina

Accurate numerical simulation of radiation reaction effects
in strong electromagnetic fields
N.V. Elkina,1 A.M. Fedotov,2 C. Herzing,1 and H. Ruhl1
1Ludwig-Maximilians Universita¨t, Mu¨nchen, 80539, Germany∗
2National Research Nuclear University MEPhI, Moscow, 115409, Russia
4
TheLandau-Lifshitzequationprovidesanefficientwaytoaccountfortheeffectsofradiationreac-
1
0 tion without acquiring the non-physical solutions typical for theLorentz-Abraham-Dirac equation.
2
We solve the Landau-Lifshitz equation in its covariant four-vector form in order to control both

n
the energy and momentum of radiating particle. Our study reveals that implicit time-symmetric
a
J collocation methods of the Runge-Kutta-Nystro¨m type are superior in both accuracy and better

0
maintaining the mass-shell condition than their explicit counterparts. We carry out an extensive
3
study of numerical accuracy by comparing the analytical and numerical solutions of the Landau-

]
h Lifshitz equation. Finally, we present the results of simulation of particles scattering by a focused
p
laserpulse. Duetoradiationreaction,particlesarelesscapableforpenetrationintothefocalregion,
-
m
ascomparedtothecaseofradiationreactionneglected. Ourresultsareimportantfordesigningthe
s
a forthcoming experimentswith high intensity laser fields.
l
p
. PACSnumbers: 41.75.Ht,41.75.Jv,03.50.De,02.60.Cb
s
c Keywords: strong laser field, classical particle dynamics, radiation reaction, Landau-Lifshitz equation, im-
i
s plicitRunge-Kutta-Nystr¨om methods
y
h
p
[

1 I. INTRODUCTION
v
1
8 In recent years much interest has been attracted to the problem of radiation reaction in a realm of high intensity
8
7 laser fields. Rigorous approach to radiation reaction requires quantum treatment of the stochastic events of hard
.
1 photon emission [1, 2]. In numerical simulation this was incorporated by using Monte-Carlo sampling of photon
0
4 emission [3]. This method treats correctly the recoil due to emission of hard photons with ~ω > mc2 0.511 MeV
∼
1
(m is the electron mass, and c is the speed of light) but is computationally rather expensive. At the same time,
:
v
i the relatively soft part of the photon spectrum 10keV - 0.511MeV had to be ignored in plasma simulations due to a
X
limited numerical resolution.
r
a
One of the approaches to account for radiation effect due to soft photon emission is to consider it classically by
addingtheradiationfrictiontermintothe equationofparticlemotion. Theclassicalproblemofmotionofaradiating
charge was considered in [4] where the so-calledLorentz-Abraham-Dirac(LAD) equation was derived. This equation
...
of motion contains third order derivative xµ and therefore possesses problematic run-away or causality violating
solutions. To ensure physical behaviour the LAD equation requires careful choice of both initial and boundary
conditions (see [5], [6] for review) which is difficult to realize in computational practice. In the most of physically
∗Electronicaddress: [email protected]
2
meaningfulsituationsradiationreactioncanbetreatedasaperturbation. Thefirstorderofperturbationtheoryresults
in the Landau-Lifshitz (LL) equation [7]. This equation contains nonlinear dissipation which leads to contraction of
a phase space volume. Proper numerical method should capture this effect.
In plasma simulation solvers [8] the equations of particle motion are usually integrated by the modified leap-frog
(St¨ormer-Verlet) method. This originally semi-implicit method can be turned to an explicit one provided that the
special ansatz [9, 10] is used to treat properly the rotation in a magnetic field. However, inclusion of an extra term
responsible for the radiation reaction breaks the symmetry of the relativistic Newton-Lorentz (NL) equation. Since
radiation reaction depends upon momenta nonlinearly, the explicitness of the method can be no more maintained
within the second order of accuracy. To avoid the loss of accuracy, one may consider different higher order methods,
primarily the single-step Runge-Kutta (RK) like methods [11].
AnotherconcernwithsimulationofthecovariantLLequationismaintainingthemassshellcondition(conservation
oftheMinkowskinormu uµ =c2). TheMinkowskinormcanbepreservedexactlyinacourseofnumericalsimulation
µ
if a corresponding equation of motion is lifted onto the Lie group settings [12, 13]. The Lie group methods show
promising results, howeverarerather expensive for nonlinear equations, as they require multiple computations of the
matrix exponentials for each time step [14].
In this paper we would like to draw attention to the Runge-Kutta-Nystr¨om [RK(N)] type methods [15], which are
speciallydesignedforsolvingthe secondorderordinarydifferentialequations(ODEs). However,the explicitRKtype
methods do not conserve exactly the quadratic integrals of motion (e.g., energy and the Minkowski norm) and hence
are not enough suitable for integration over a long time. Contrarily, the implicit collocation RK(N) methods [16]
provide better accuracy and are capable for conservation of quadratic integrals of motion. This, together with time
reversibility, singles out these methods as perfect candidates for integration of the LL equation over a long time.
The paper is organized as follows. In Section II a summary of covariant formulation of the equation of motion of
radiating particle is presented. The main goal of the paper is achieved in Section III, where we present derivation
of the collocation Gauss-Legendre RKN method entirely in terms of collocation points for fourth, sixth and eighth
orders of accuracy. The result is employed in Sections IV and V, where we perform extensive validation studies of
the consideredRK(N)methods againstexactsolutionsofthe LL equationin aconstantmagneticfield andin aplane
wave and demonstrate that implicit collocation methods of RK(N) type are superior in accuracy to their explicit
counterparts in solution of both the NL and LL equations. In Section VI we present an important application of
the developed numerical technique to the problem of relativistic scattering of electrons by a focused laser field. Our
conclusions are summarized in Section VII.
II. BASIC EQUATIONS
Radiation reaction modifies the NL equation by a radiation friction term gµ in the right hand side
d2xµ e
=fµ = Fµuν +gµ. (1)
dτ2 m ν
3
The 4-velocity uµ = γ, u of a charge is the derivative of the position 4-vector xµ = t, r with respect to the
{ } { }
proper time τ. The electromagnetic field tensor Fµ is expressed by the following matrix
ν
0 E E E
x y z
 
E 0 B B
Fµ = x z − y 
ν  
 E B 0 B 
 y − z x 
 
 
 E B B 0 
 z y x 
 − 
 
Radiation reaction was originally considered with the LAD equation [4]
d2xµ e 2e2 d2uµ d2u
= Fµuν + u uµ ν uν, (2)
dτ2 m ν 3m ν dτ2 − dτ2
(cid:18) (cid:19)
here and in the following we adopt the units c = ~ = 1, e = √α, α = e2/~c 1/137 is the fine structure constant.
≃
As was pointed out in the Introduction, the equation (2) possesses the problematic run-away or causality violating
solutions. The regular way to avoid such unwanted solution is to treat the second term in right-hand-side as a
perturbation [7]. Then in virtue of
dxλ∂Fµ
F˙µ = ν =uλ∂ Fµ
ν dτ ∂xλ λ ν
the term gµ takes the form
2e3 2e4
gµ = uλ(∂ Fµ)uν + FµFνuλ+(F ul)(Fkum)uµ . (3)
3m2 λ ν 3m3 ν λ kl m
(cid:2) (cid:3)
Eq.(1)withthe termgµ givenby(3)isknownasthe LLequation. The numericalsolutionofthe LLequationshould
preserve the Minkowski norm uµu = 1, as well as the orthogonality condition gµu = 0. It turns out that actual
µ µ
fulfilmentoftheformerconditionstronglydepends onachosennumericalmethod,while implementationofthe latter
one is difficult to enforce if gµ is in the form (3). Indeed, in this form gµ consists of three terms, which contribute
disparately. While in the majorityof situations the firstterm in (3) can be safely neglected, the two remainingterms
together do satisfy the orthogonality condition, but due to a huge difference in the magnitude of their components
naive summation would lead to accumulation of considerable round-off errors. This purely numerical problem may
cause eventual violation of orthogonality between the 4-force and 4-momentum. In order to avoid such an unwanted
effect we consider the LL equation in the “symmetrized” form reminiscent of the LAD equation. So that the actual
equation which is solved in this paper is given as follows
d2xµ e 2e2
= Fµuν + (u w˙µ uµw˙ )uν, (4)
dτ2 m ν 3m ν − ν
where w˙µ is calculated as
e e e
w˙µ = F˙µuν +Fµu˙ν = uλ(∂ Fµ)uν + FµFνuλ .
m ν ν m λ ν m ν λ
(cid:16) (cid:17) h i
InEq.(4)both leadingradiationreactionterms acquirethe sameorder,sothattheir subtractioncannolongercause
accumulation of round-off errors. Such a trick allows to keep gµu =0 satisfied up to the level of machine precision.
µ
4
III. NUMERICAL INTEGRATION OF THE MOTION OF A CHARGE
A. Explicit Runge-Kutta-Nystro¨m methods
The mostfrequently preferredapproachto solvea secondorderODE x¨µ =fµ(τ,xµ,x˙µ) with higher accuracyis to
apply a suitable method to an equivalent system of the first order system of ODEs
d xµ uµ
 = . (5)
dτ
uµ fµ(τ, xµ, uµ)
   
   
   
Application of a RK method to such a problem results in [25]
s
Lµ =uµ+h A Kµ, i=1,...,s, (6)
i n ij j
j=1
X
s s
Kµ =fµ τ +hc , xµ+h A L , uµ+h A K , i=1,...,s, (7)
i  n i n ij j n ij j
j=0 j=1
X X
 s 
xµ =xµ+h a Lµ, (8)
n+1 n j j
j=1
X
s
uµ =uµ+h a Kµ, (9)
n+1 n j j
j=1
X
where h is a step size and A , c and a are the scheme-specific coefficients (see Table Ia). An alternative approach
ij i i
TABLE I: The Butcher tables for RK and RKN methods.
c A c A B
s s s s s s
× × ×
aT aT bT
(a)RK (b)RKN
tosolvethe equationofmotionis toexplorethe originalsecondorderODEx¨µ =fµ(τ,xµ,x˙µ)withoutpriorreducing
it to a system of the first order ODEs. This can be achieved by eliminating L from the formulas (6) - (9) as follows
i
s s
Kµ =fµ τ +c h,xµ+c huµ+h2 B Kµ,uµ+ A Kµ , (10)
i  n i n i n ij j n ij j
j=1 j=1
X X
 s 
xµ =xµ+huµ+h2 b Kµ, (11)
n+1 n n j j
j=1
X
s
uµ =uµ+h a Kµ, (12)
n+1 n j j
j=1
X
where the coefficients B and b are given by
ij j
B = A A , b = a A . (13)
ij ik kj j k kj
k k
X X
5
The conditions (13) for the new coefficients B , b in fact can be abandoned, as it was found by Nystr¨om in [15].
ij j
Such kind of methods is known as the RKN methods, see also [17]. In most general form the RKN methods are
defined by Eqs. (10) - (12), where the coefficients A , B , a , b are kept independent and form the double Butcher
ij ij i i
array (see Table (Ib)). If A = 0 and B = 0 for j i, then such construction yields an explicit method, in which
ij ij
≥
TABLE II: The explicit RK(N) methods of order 4.
0 0 0 0 0 0 0 0 0 0 0 0 0 0
1/2 1/2 0 0 0 1/2 1/2 0 0 0 1/8 0 0 0
1/2 0 1/2 0 0 1/2 0 1/2 0 0 1/8 0 0 0
1 0 0 1 0 1 0 0 1 0 0 0 1/2 0
1/6 1/3 1/3 1/6 1/6 1/3 1/3 1/6 1/6 1/6 1/6 0
(a)RK (b)RKN
eachstage Kµ can be determined sequentially. The examples of explicit forth order RK and RKN methods are given
i
in Table II.
B. Implicit collocation Runge-Kutta-Nystro¨m methods
The basic idea of the collocation methods [16] is to select such polynomials rµ(τ) of degree s+1 that satisfy the
initial value problem at s+1 points τ , τ +c h, ..., τ +c h [τ ,τ ), where c ,...,c are the distinct real
n n 1 n s n n+1 1 s
∈
numbers from the interval (0, 1). Then for the second order ODE x¨µ = fµ(τ,xµ,x˙µ) the corresponding polynomial
approximation rµ satisfies [18]
rµ(τ )=xµ, r˙µ(τ )=uµ, r¨µ(τ +c h)=fµ(τ +c h, rµ(τ +c h), r˙µ(τ +c h)) [i=1,...,s]. (14)
n n n n n i n i n i n i
In the following we proceed with a derivation of an RKN method (the corresponding RK method is defined by
a subtable with the coefficients A ). Being of degree s 1, the second derivative of the collocation polynomial
ij
−
r¨µ(τ +ξh) can be uniquely expanded in the Lagrange polynomials
n
s
ξ c
k
L (ξ)= − . (15)
i
c c
i k
k6=i −
Y
Setting Kµ =r¨µ(τ +c h), this expansion can be written as
i n i
s
r¨µ(τ +ξh)= L (ξ)Kµ, τ +ξh [τ , τ +h].
n j j n ∈ n n
j=1
X
This implies that r˙µ and rµ can be approximated on [τ , τ +h] by
n n
s ξ
r˙µ(τ +ξh)=uµ+ KµA (ξ), A (ξ)= L (ξ′)dξ′, (16)
n n i i i j
i=1 Z0
X
s ξ ξ
rµ(τ +ξh)=xµ+ξh uµ+ KµB (ξ), B (ξ)= A (ξ′)dξ′ = (ξ ξ′)L (ξ′)dξ′. (17)
n n · n i i i i − i
i=1 Z0 Z0
X
6
Now,bysubstitutingξ =c forsuccessivei=1,...,s,wecanmaketheformulas(14),(16),(17)coincidingtoEqs.(10)
i
- (12), and thus specify the discrete coefficients of the collocation method:
A =A (c ), a =A (1), B =B (c ), b =B (1). (18)
ij j i j j ij j i j j
The forth order method can be obtained by using two Lagrange polynomials of degree one
ξ c ξ c
2 1
L (ξ)= − , L (ξ)= − , (19)
1 2
c c c c
1 2 2 1
− −
the resulting method is presented in Table IIIa. The highest possible order 2s can be achieved when the collocation
TABLE III: Double Butcher table for the collocation RKN method for s=2.
c1 c12((cc11−2c2c2)) 2(c1c21c2) c216((cc11−3c2c2)) 3(c1c31c2) 21 − √63 14 14 − √63 316 356 − √123
− − − −
c2 2(c2c22c1) c22((cc22−c21c1)) 3(c2c32c1) c226((cc22−3c1c1)) 12 + √63 14 + √63 41 356 + √123 316
− − − −
2(1c−12cc22) 2(1c−22cc11) 6(1c−13cc22) 6(1c−23cc11) 12 21 14 + √123 41 − √123
− − − −
(a)CollocationRKN (b)Gauss-LegendreRKN
points are chosen as roots of the sth shifted Legendre polynomial [11]
ds
[τs(τ 1)s]=0, τ (0, 1). (20)
dτs − ∈
Other possible choices of the collocation points (Gauss-Lobatto, Gauss-Radau methods) possess lower order p < 2s,
but might be still beneficial if e.g. the equations are stiff [19]. For s = 2 the roots of the quadratic polynomial (20)
are c =1/2 √3/6 and the resulting Gauss-Legendre-RKNmethod of order 4 is presented in Table IIIb.
1,2
∓
In this paper we also use higher order methods based on s = 3 and s = 4 collocation points. The Gauss-
Legendre-RKN method with s = 3 collocation points is briefly described in Appendix A. Though computationally
more expensive, these 6th and 8th ordermethods provide possibility to measurethe accuracyoflower ordermethods
in case analytical solution is not available.
After thecoefficients(18)aresubstitutedintothe system(10)-(12)itturns toasetofstageequations. Implicitness
of the collocationRK(N) methods implies usage of the iterationtechnique to solve the stage equations. The simplest
iteration scheme can be obtained by just hanging the iteration induces m, m+1 on the implicit equations:
s s
iRK(N): Kµ,m+1 =fµ(τ +c h, xµ+c huµ+h2 B Kµ,m, uµ+h A Kµ,m), (21)
i n i n i n ij j n ij j
j=1 j=1
X X
where the prefix “i”stands for “implicit” to distinguishthe method fromexplicitones, for the latter we will similarly
use the prefix “e” in front of their abbreviation RK(N). In order to formulate the criterion for stopping the iteration
process, an error of the approximation has to be measured by an appropriate norm. In this paper we use for this
purpose the maximum norm
Km+1 Km =max Kµ,m+1 Kµ,m . (22)
|| − ||∞ i,µ | i − i |
In the rest of the paper, we discuss applications of the described numerical methods to physical problems of particle
motion in ultrastrong electromagnetic fields.
7
IV. PARTICLE IN MAGNETIC FIELD
Let us start with motion of an ultrarelativistic particle in a constant magnetic field. By fixing the direction of
magnetic field along x3 =z, the LL equation reads
du e 2 e4
⊥ = u H H2(1+u2)u , (23)
dτ m ⊥× − 3m3 ⊥ ⊥
du 2 e4
z = H2u2u , (24)
dτ −3m3 ⊥ z
where u = u , u is the transverse component of 4-velocity. After passing to polar coordinates u = u cosφ,
⊥ x y x ⊥
{ }
u = u sinφ in the transverse plane xy and introduction of dimensionless time ϕ = Ωτ (where Ω = eH/m is the
y ⊥
rotation frequency in the proper reference frame), the system (23) - (24) takes the form
du
⊥ = Ku (1+u2) (25)
dϕ − ⊥ ⊥
dφ
= σ, (26)
dϕ −
du
z = Ku2u . (27)
dϕ − ⊥ z
Here σ = 1 is the sign of a particle charge and K is the dimensionless parameter
±
2e2Ω 2e3H 2αH
K = = = , (28)
3 m 3 m2 3 H
c
whereH =m2/eis the criticalmagnetic field. Note thatunder applicabilityofLL equationK 1. Eq.(25)admits
c
≪
separation of variables, after that the rest equations (26) - (27) can be also easily solved. The solution reads
1
u(ϕ)= u cos(φ σϕ), u sin(φ σϕ), u eKϕ , (29)
⊥0 0 ⊥0 0 z0
(1+u2 )e2Kϕ u2 − −
⊥0 − ⊥0
(cid:8) (cid:9)
wherethe subscript“0”refperstotheinitialvaluesatϕ=0. Bysimpleintegrationonecanfindtheparticletrajectory
in terms of hypergeometric function, e.g. its projection onto the plane xy is given by
ϕ
x x0 σ 1 iσ iσ (1+u2 )
 y−−y0 = u⊥0 q(1+u2⊥0)e2Kϕ−u2⊥0−ℑℜei(φ0−σϕ)·2F1(cid:18)1, 2 − 2K, 1− 2K, u2⊥0⊥0 e2Kϕ(cid:19)(cid:12)(cid:12)(cid:12)0,
(cid:12)
     (cid:12)(30)
(below we assume u = 0). The trajectory is a converging spiral shown in Fig. 1a. Trajectory without radiation
z0
reaction corresponds to the case K = 0 and represents a circle. Numerical integration is performed with the RK
method, since this particular problem does not contain dependence on coordinates. The observed order of accuracy
of the method can be estimated by examining the behaviour of L norm
2
h
L = uµ(τ ) uµ (τ )2, (31)
2 i ex i
T | − |
s i
X
where T is the proper time duration of simulation and uµ refers to the exact solution (29).
ex
The results of studies of convergence of different numerical methods are presented in Fig. 1b, where both L and
2
the Minkowskinormerrorsareplottedindependenceuponatime steph. Theslopesofthecurvesintheupperpanel
8
10−1 eRK4
iRK4
iRK6
2 10−5 iRK8
L
4 6
∝ h ∝ h 8
10−9 h
∝
10−1
|
1
− 10−5
µ
u
µ 10−9
u
|
10−13
10−5 10−4 10−3 10−2
time step, h
(a)Particletrajectories. (b)Ontop: dependence oftheerrormeasureL2 uponatimestep;
onbottom: inaccuracyinpreservationoftheon-shellcondition.
FIG. 1: (Color online). Effect of the radiation reaction on particle motion in a constant magnetic field
eH /mΩ=103. Initial γ-factor of a particle is γ =103.
z 0
of Fig. 1b, corresponding to both the explicit eRK4 and implicit iRK4 methods clearly correspond to the order four,
as expected, but the difference in the initial error level makes the implicit method more favourable with respect to
higher accuracy. The higher order methods iRK6 and iRK8 also exhibit the expected slopes, corresponding to orders
six and eight, respectively. The inaccuracy of preservation of the mass shell condition, as shown in the lower panel
of Fig. 1b, remains low and bounded for all implicit methods, while the eRK method results in substantial violation
of the Minkowski norm. It should be noted that non-conservation of the Minkowski norm depends on the numerical
error δuµ as
(uµ+δuµ)(u +δu )=1+2uµδu +O(δu2),
µ µ µ
and hence is especially pronounced in ultrarelativistic case uµ 1.
≫
V. PARTICLE IN PLANE WAVE FIELD
The LL equation admits much simpler form for a particle moving in a plane wave field A (x)=A (ϕ), ϕ=kx=
µ µ
ω t k r, kA=0:
0
− ·
duµ e 2 e3 2 e4 2 e4
= h(A′u)kµ Aµ′ + h(A′′u)kµ Aµ′′ A′2kµ+ A′2uµ (32)
dϕ m − 3m2h − − 3m3 3m3h
(cid:2) (cid:3) (cid:2) (cid:3)
9
Newton-Lorentzequation Newton-Lorentzequation
Landau-Lifshitzequation Landau-Lifshitzequation
2500 2500
2000 2000
1500 1500
0 0
ω ω
/ 1000 / 1000
c c
500 500
, ,
y y
0 0
−500 −500
−1 0 1 2 3 4 5 −1 0 1 2 3 4 5
˜ ˜
x, c/ω x, c/ω
0 0
(a)Circularpolarization (b)Linearpolarisation
FIG. 2: (Color online). Trajectoriesof a chargein a plane wavefield for a =103 in a frame that wouldbe comoving
0
for a guiding center in the absence of radiation reaction. Radiation reaction leads to a drift of guiding center along
the wave vector.
where1/h=ϕ˙ =k uµ. Thisequationcanbesolvedanalyticallyinquadratures[20]andinspecialcasesthe4-velocity
µ
uµ can be even expressed in terms of elementary functions [21]. Here we consider two important cases corresponding
to linear and circular polarisation of the wave.
The effect of radiation reaction on a particle trajectory is particularly transparent in a reference frame comoving
with the guiding center, as in this frame the radiation-free trajectories take especially simple form: in case of linear
polarization particles move along the figure-8, whereas in case of circular polarization the trajectory is a circle. The
radiationreactioneffects openclosedfigure-8orcircular motionandleadto a drift alongthe wavevectorofthe wave
(see Fig. 2).
Inordertogainmoreknowledgeabouttheaccuracyofthenumericalmethodswecomparenumericalandanalytical
solutions in case of a circularly polarized plane wave, for which the dependence of the phase ϕ and therefore of uµ
and xµ on the proper time τ can be expressed explicitly by [21]
1+2a2(ku )2Kτ/ω 1
φ(τ)= 0 0 0− , (33)
a2(ku )K/ω
p 0 0 0
where K is this time defined by K =(2/3)e2ω /m. The results are presented in Fig. 3. One can see that the slopes
0
of the curves in Fig. 3a are in agreement with the expected orders of the methods, indicated by dashed lines. Only
the curveforeRKN4 method onthis plotdeviates slightlyfrompowerlaw whenthe time stepbecomes larger. As for
Fig. 3b, it demonstrates the same behavior of the methods as was previously pointed out for the problem of motion
in a constant magnetic field (compare to Fig. 1b).
Let us discuss computationalcomplexity of the implicit methods. The mainsource ofcomputationalcomplexity in
implicit methods comes from the number of iterations needed for a nonlinear solver to achieve the desired accuracy.
Hence we measure it for both the RK and RKN methods by an average number of iterations per step. The results
for a case of motion in the field of linearly polarized plane wave are summarized in Fig. 4.
According to Fig. 4, the iRKN method requires less iterations than iRK. One can also compare qualitatively the
10
eRKN4 eRKN4
iRKN4 iRKN4
10−1 iRKN6 iRKN6
iRKN8 101 iRKN8
h2
∝
m 1| 10−2
r −
o
n µ
u
L2 10−5 µ
u
|
h4 10−5
∝
h8
∝
h6
∝ 10−8
10−9
10−2 10−1 100 10−2 10−1 100
time step, h time step, h
(a)Convergence inL2 norm (b)ErrorinMinkowskinorm
FIG. 3: (Color online). Convergence study for a particle moving in a circularly polarised plane wave for a =103.
0
performance of implicit iRKN4 and its explicit counterpart eRKN4. The number of stage equations (10) for eRKN4
equals 4,while the correspondingnumber ofstages for iRKN4 is 2. But because of4 iterations onthe averagefor the
same time step, the iRKN4 method is twice more expensive in the number of calls of the right-hand-side function.
From the other hand the implicit methods admit longer time steps without deterioration of numerical accuracy.
Moreover,furtherincreaseoftheorderofexplicitmethods eventuallyresultsinreachingtheso-calledButcherbarrier
[11], after which the amount of stages starts to grow faster than the order.
VI. SCATTERING OF PARTICLE ON FOCUSED LASER PULSES
In this section we apply the developed numerical technique to an important problem of scattering of particles by
focusedlaserpulses. Intheabsenceofradiationreaction,thisproblemcanbestudiedwiththetheoryofponderomotive
scattering, see e.g. [22, 23]. In previous studies of this problem radiation reaction was neglected because of the
assumption of relatively small both the intensity and the particle energy. Our basic settings are similar to those
considered in [23] but with the laser intensity and particle energy higher by two order of magnitude, γ , a 103.
0 0
∼
Consider first a focused monochromatic Gaussian beam with the electric field defined by
iE b2 ω[r (ε n)](r n) ω(r n)2
E = 0 ε+ · × × e−iω0(t−r·n)exp × . (34)
(b+ir n)2 b+ir n −2(b+ir n)
· (cid:26) · (cid:27) (cid:20) · (cid:21)
Our notations here are the same as in Ref. [24]: ω is the carrier frequency of the beam, b is a half of the Rayleigh
0
length, ε is a (complex) polarization vector, and n is a unit vector along the focal axis. The magnetic field is given
byH = (ic/ω ) E. Ourtargetconfigurationconsistsoftwobeamswiththe samepolarizationε,but oppositely
0
− ∇×
directed vectors n. In such a case near the focal axis the electric field is doubled, while the magnetic field vanishes.