J Dubois

S A Berman

C Chandre

T Uzer

Guiding-center motion for electrons in strong laser fields

We consider the dynamics of electrons in combined strong laser and Coulomb fields. Under a time-scale separation condition, we reduce this dynamics to a guiding-center framework. More precisely, we derive a hierarchy of models for the guiding-center dynamics based on averaging over the fast motion of the electron using Lie transforms. The reduced models we obtain describe well the different ionization channels, in particular, the conditions under which an electron is rescattered by the ionic core or is directly ionized. The comparison between these models highlights the models which are best suited for a qualitative and quantitative agreement with the parent dynamics.

I. INTRODUCTION

Electrons in atoms or molecules subjected to laser fields can tunnel ionize when the laser intensity is strong enough to compete with the Coulomb attraction [1][2][3][4]. After tunneling, the electron can for example ionize directly, or be rescattered by the ionic core [4][5][6]. Rescattering occurs when the electron tunnel ionizes, then comes back to the ionic core region, and is finally ionized after a strong interaction with the ionic core. This process, also referred to as a recollision, is the keystone of strong field physics [7]. In contrast, direct ionization occurs when the electron never comes back to the ionic core after tunnel ionization. The parameters of the laser pulse, e.g., its intensity, frequency, polarization, and envelope, play a paramount role in the ionization process, as reflected by the shaping in the momentum distributions of the ionized electrons measured at the detectors [8][9][10][11]. In order to better analyze and interpret what is measured at the detectors, there is a need to better understand the motion of the electrons after tunneling [12][13][14][15][16][17].

To describe the motion of the electron after tunneling, a convenient approximation is usually performed: Since the electron is relatively far from the ionic core after tunneling and the electric field is strong, the Coulomb interaction is neglected. This leads to the so-called strong field approximation [5] (SFA), which is one of the main theoretical tools in strong-field atomic physics. The main advantage of this approximation is that it allows the explicit computation of the trajectories of the electron since the equations of motion are linear. However, the SFA often leads to disagreements or misleading interpretations when it is confronted with experimental data [10,[18][19][20]. For example, in linearly polarized (LP) fields, the SFA suggests that if the electron does not return to the ionic core within one laser cycle after tunneling, it never comes back to the ionic core at all. However, multiple laser cycle recollisions are essential for the quantitative agreement between theories and experiments in non-sequential double ionization (NSDI) [18,[21][22][23][24]. In circularly polarized (CP) fields, the SFA suggests that the drift-velocity of the electron pushes it away from the ionic core, without rescattering. The predicted absence of recollisions is in contradiction with the knee structure observed in the double ionization probability curves, a signature of the recollision process, observed in simulations [25][26][27] and experiments [20]. In fact, multiple laser cycle recollisions or recollisions in CP become possible only when the Coulomb interaction is taken into account. In addition, Coulomb effects such as Coulomb focusing [18,28] and Coulomb asymmetry [19,29] are measurable and significant in above-threshold ionization (ATI) [10,28,30]. Therefore, the Coulomb potential cannot be ignored, even far away from the ionic core.

In the length gauge (see Ref. [31] for a review), the Hamiltonian of an electron interacting with its parent ion and a laser field reads

H(r, p, t) = |p| 2 2 + V (r) + r • E(t), (1)
where r and p are the position of the electron and its canonically conjugate momentum, respectively. Atomic units (a.u.) are used unless stated otherwise. Here, we have considered a single-active electron atom, where the effective charge is equal to one due to the shielding effect by the inner electrons. Furthermore, the timescale of interest is of the order of femtoseconds. During this short timescale, the ionic core is too massive to display relevant variations of its position, and consequently, we have considered the ionic core to be static. Moreover, the characteristic distance of the electron in the laser field is small compared to the wavelength of the laser field, and consequently, we have neglected the spatial dependence of the electric field (dipole approximation). The polarization of the electric field is crucial in the analysis of the variety of nonlinear phenomena [26,[32][33][34][35][36][37][38][39][40]. The elliptically polarized electric field is

E(t) = f (t) E 0 ξ 2 + 1 [x cos(ωt + φ) + ŷξ sin(ωt + φ)] .
(2) The laser ellipticity is ξ, where in LP and CP fields, ξ = 0 and ξ = 1, respectively. The laser amplitude in a.u. is given by E 0 ≈ 5.338 × 10 -9 √ I, with I the laser intensity in W • cm -2 . The laser phase at t = 0 is φ. In this manuscript, we consider the laser envelope to be f (t) = 1. Here, we consider a single frequency ω for the laser field.

Figure 1 displays six typical trajectories of Hamiltonian (1). We notice that these trajectories display fast oscillations around a mean trajectory, which we call the guiding-center trajectory. In this article, we use this clear separation of scales to derive models for the guidingcenter dynamics.

Specifically, we average Hamiltonian (1) over the fast motion using Lie transforms in order to simplify the electron dynamics and clearly distinguish the different ionization channels for the electron. As a consequence, we derive a hierarchy of averaged models which fully take into account the Coulomb potential. In these models, the electron in the combined strong laser and Coulomb fields follows a guiding center trajectory [30]. Actually, there are several possible guiding centers, depending on the order up to which the averaging is performed. Computations are presented for a d-dimensional configura-tion space, where d = 1 (only for linear polarization) or d = 2, 3. The objective of this manuscript is to derive and investigate these reduced models.

In Sec. II, we describe the procedure we use for averaging Hamiltonian (1) over the fast time scale, and we derive a hierarchy of models for the guiding center dynamics. In Sec. III, we first compare these reduced models with the dynamics associated with Hamiltonian (1). In particular, we show the relevance of two models in the hierarchy. Then we analyze the dynamics of the guidingcenter models in phase space, highlighting regular and chaotic regions and their relation with the trajectories of Hamiltonian (1).

II. DERIVATION OF THE HIERARCHY OF MODELS A. Reminders on canonical Lie transforms

In this section, we recall some basic features on canonical transformations in the framework of canonical Lie transforms. For more details, we refer to Refs. [41,42]. We consider a Hamiltonian system with phase-space variables z, a Hamiltonian H(z) and a Poisson bracket {•, •} (whether it is a canonical or non-canonical bracket). Canonical Lie transforms are near-identity changes of coordinates z → z(z), generated by a scalar function S(z), called the generating function, and given by z = exp (-L S ) z, for any scalar functions F (z) and G(z). As a consequence of Property (i) and of the scalar invariance F (z) = F (z), these changes of variables modify any observable F (z), and in particular the Hamiltonian H(z), into

F (z) = exp (L S ) F (z). (4)
Property (ii) ensures that these changes of coordinates do not affect the expression of the Poisson bracket, i.e., they are canonical transformations. One significant advantage of these transformations is that they are explicit functions and they can be easily inverted: z = exp(L S )z. This way, we can recover all the information on the particle dynamics from the transformed (averaged) system. These canonical Lie transforms are particularly well suited for perturbation theory. If the Hamiltonian is of the form H = H 0 + ǫW , where H 0 is the Hamiltonian of the unperturbed system, W is the perturbation and ǫ is an ordering (small) parameter, a canonical Lie transform generated by a generating function S(z) (of order ǫ), applied to H is able to remove the unwanted part of the perturbation, called W , and move its influence to higher orders in ǫ. More explicitly, at the lowest order, the expression of the new Hamiltonian expressed in the new variables z is

H = exp (L S ) H (5) = H 0 + ǫW + {S, H 0 } +ǫ{S, W } + 1 2 {S, {S, H 0 }} + . . . .
Choosing appropriately the generating function S such that {S, H 0 } = -ǫ W , unwanted terms in the perturbation W can be pushed from order ǫ to order ǫ 2 , meaning that the order ǫ in the Hamiltonian becomes W -W . For instance, one can suppress fast oscillating terms contained in W . Then, the associated canonical change of coordinates is determined using Eq. (3).

B. Averaging the electron dynamics

We notice that these transformations are defined for autonomous systems. Hamiltonian (1) has an explicit time dependence through the electric field. Therefore, we first increase phase space to include time t, and consider its canonically conjugate variable k. The extended Hamiltonian (1) becomes

H(r, p, t, k) = k + |p| 2 2 + V (r) + r • E(t). (6)
The extended Poisson bracket is

{F, G} = ∂F ∂r • ∂G ∂p - ∂F ∂p • ∂G ∂r + ∂F ∂t ∂G ∂k - ∂F ∂k ∂G ∂t , (7)
where the operators ∂/∂r = (∂/∂x, ∂/∂y, ∂/∂z) and ∂/∂p = (∂/∂p x , ∂/∂p y , ∂/∂p z ). The hypothesis we make for the derivation of our hierarchy of reduced models is that the characteristic time of the ionized electron trajectory is large compared to a laser cycle T = 2π/ω, i.e., ω → ω/ǫ where ǫ is an ordering parameter which is explicitly introduced for bookkeeping purposes. Performing the canonical change of coordinates t = t/ǫ and k = ǫk, and re-scaling the energy, Hamiltonian (6) becomes

H (0) (r, p, t, k; ǫ) = k + ǫ |p| 2 2 + V (r) + r • E(t) , (8)
where we have removed the bars in the new variables. We apply canonical Lie transforms as described above in order to perform the averaging of Hamiltonian (8) over the fast time scale, by pushing time-dependent terms in the Hamiltonian to higher order terms in ǫ.

Gauge velocity transformation

As an example, we consider the transformation from the length gauge to the velocity gauge [31] in Hamiltonian (8). It is given by the following change of coordinates r = r, p = p -ǫA(t), where A(t) is the vector potential defined by E(t) = -∂A(t)/∂t. This transformation is a canonical change of coordinates which can be formulated as a canonical Lie transform generated by

S (1) = ǫ r • A(t). (9)
The Hamiltonian in the velocity-gauge coordinates becomes

H (1) = exp(L S (1) )H (0) (10) = k + ǫ 1 2 (p + ǫA(t)) 2 + V (r) , = k + ǫ |p| 2 2 + V (r) + ǫ 2 p • A(t) + ǫ 3 A 2 (t) 2 .
We observe that the time-dependence in Hamiltonian ( 8), present in E(t), is of order ǫ, while in Hamiltonian (10) this time-dependence is moved to order ǫ 2 .

Iterative procedure

We iterate the above-procedure to higher order in ǫ. We assume that after the N -th step of the procedure, all the time-dependent terms in the present averaged Hamiltonian are removed up to order ǫ N , that is, the time-dependence of the averaged Hamiltonian is of order ǫ N +1 . We assume that the total generating function up to order ǫ N is known. The total generating function at this step is

S (N ) (r, p, t; ǫ) = N n=1 ǫ n S n (r, p, t), (11)
and the corresponding averaged Hamiltonian is denoted

H (N ) = exp(L S (N ) )H (0) = k + N n=1 ǫ n h n (r, p) + ǫ N +1 R N +1 (r, p, t; ǫ),
where h n (r, p) are the coefficients in the series expansion of the Hamiltonian that no longer depends on time, while R N +1 (r, p, t; ǫ) is the remainder of the Hamiltonian which still depends on time. The objective of the iterative method is to find the modified generating function S (N +1) [which amounts to finding the extra function S N +1 in Eq. (11)] to remove the time-dependence in the term R N +1 at the lowest order. The averaged Hamiltonian H (N +1) whose timedependence is of order ǫ N +2 is

H (N +1) = exp(L S (N +1) )H (0) , = exp(ǫ N +1 L SN+1 ) exp (L S (N ) ) H (0) + O(ǫ N +2 ), = exp(ǫ N +1 L SN+1 )H (N ) + O(ǫ N +2 ), = H (N ) + ǫ N +1 R N +1 + ∂S N +1 ∂t + O(ǫ N +2 ).
The time-fluctuating terms in R N +1 are denoted R N +1 , and are defined by

R N +1 = R N +1 - 1 T T 0 dtR N +1 .
In order to eliminate the time-fluctuating terms at order ǫ N +1 , the component S N +1 of the generating function S (N +1) is chosen as

S N +1 = -dt R N +1 (r, p, t; 0), ( 12
)
where the primitive is chosen such that the mean value of S N +1 with respect to t ∈ [0, T ] is zero. At each step the functions R n have to be computed up to order M where M is the last order for which the averaged Hamiltonian will be computed analytically. We perform these computations using a symbolic computation software.

Averaged Hamiltonians

We apply the above-described procedure to push the time-dependence in Hamiltonian (8) to order ǫ 8 using H (7) = exp(L S (7) )H (0) . Below we provide the explicit expression for S (6) . The higher-order components are too lengthy to report and their expressions are not particularly enlightening.

S (6) = ǫ r • A(t) - ǫ 2 ω 2 p • E(t) + ǫ 3 ω 2 A(t) • E(t) 4 + ∂ ∂r V - ǫ 4 ω 4 p • ∂ ∂r E(t) • ∂ ∂r V + ǫ 5 ω 4 E(t) 4 + ∂V ∂r • ∂ ∂r -p • ∂ ∂r 2 A(t) • ∂ ∂r V - 5ǫ 6 8ω 6 p • ∂ ∂r E(t) • ∂ ∂r 2 V - ǫ 6 ω 6 p • ∂ ∂r ∂V ∂r • ∂ ∂r -p • ∂ ∂r 3 + 2 ∂V ∂r • ∂ ∂r p • ∂ ∂r E(t) • ∂ ∂r V. (13)
Here, we have used the fact that the electric field is monochromatic and satisfies ω 2 E(t) = -∂ 2 E/∂t 2 . The averaged Hamiltonian H (7) is

H (7) = k + ǫ |p| 2 2 + V (r) + ǫ 3 U p + ǫ 5 U p ω 2 (ξ 2 + 1) ∂ 2 V ∂x 2 + ξ 2 ∂ 2 V ∂y 2 +ǫ 7 U p ω 4 (ξ 2 + 1) ∂ ∂r ∂V ∂x 2 + ξ 2 ∂ ∂r ∂V ∂y 2 + O(ǫ 8 ), (14)
where all the derivatives are evaluated at r. By truncating the Hamiltonian at a given order, we notice that the reduced (time-independent) Hamiltonians up to order ǫ 7 are of the form

H(r, p) = |p| 2 2 + V eff (r),
with an effective potential V eff . In particular, this highlights a particular property in the reduction process that the reduction procedure does not generate p-dependent terms in the Hamiltonian other than the kinetic energy, up to order ǫ 7 . At order ǫ 8 , the term which is generated in H (8) is linear in the momenta; therefore, it can easily be eliminated by a translation in p (which is a canonical transformation). At order ǫ 9 , the terms which are generated are quadratic in the momenta p, and we do not consider these terms in what follows. Therefore, our analysis is valid up to order ǫ 8 . In such a way, it is particularly convenient to define effective potentials. Depending on the order of truncation, we obtain three reduced Hamiltonians

H 2 (r, p) = |p| 2 2 + V (r), (15a)
H 5 (r, p) = |p| 2 2 + V (r) + U p ω 2 (ξ 2 + 1) ∂ 2 V ∂x 2 + ξ 2 ∂ 2 V ∂y 2 , (15b)
H 7 (r, p) = |p| 2 2 + V (r) + U p ω 2 (ξ 2 + 1)
∂ 2 V ∂x 2 + ξ 2 ∂ 2 V ∂y 2 + U p ω 4 (ξ 2 + 1) ∂ ∂r ∂V ∂x 2 + ξ 2 ∂ ∂r ∂V ∂y 2 , (15c)
where we have removed the small parameter ǫ which was originally introduced for bookkeeping purposes. Each of these Hamiltonians describes the dynamics of the guiding center at a different level of approximation. As a result of averaging, the Hamiltonians (15) are conserved, in contrast to Hamiltonian (1). We notice that the quantities of the effective potentials depend on the main parameters of the electric field, its amplitude and its ellipticity, with the exception of H 2 . The corresponding change of coordinates which maps H (0) (r, p, t, k; ǫ) to H (7) is given by Eq. ( 3) and its series expansion up to order ǫ 6 is given by

r = r - ǫ 2 ω 2 E(t) - ǫ 4 ω 4 ∂ ∂r E(t) • ∂ ∂r V - 2ǫ 5 ω 4 ∂ ∂r p • ∂ ∂r A(t) • ∂ ∂r V - ǫ 6 8ω 6 ∂ ∂r E(t) • ∂ ∂r 2 V - ǫ 6 ω 6 ∂ ∂r ∂V ∂r • ∂ ∂r -3 p • ∂ ∂r 2 E(t) • ∂ ∂r V - 2ǫ 6 ω 6 ∂V ∂r • ∂ ∂r E(t) • ∂ ∂r ∂V ∂r + O(ǫ 7 ), (16a)
p = p -ǫA(t) - ǫ 3 ω 2 ∂ ∂r A(t) • ∂ ∂r V + ǫ 4 ω 4 ∂ ∂r p • ∂ ∂r E(t) • ∂ ∂r V - ǫ 5 ω 4
∂ ∂r

E(t) 4 + ∂V ∂r • ∂ ∂r -p • ∂ ∂r 2 A(t) • ∂ ∂r V + ǫ 6 8ω 6 p • ∂ ∂r E(t) • ∂ ∂r 2 ∂V ∂r + ǫ 6 ω 6 ∂ ∂r p • ∂ ∂r ∂V ∂r • ∂ ∂r -p • ∂ ∂r 3 + 2 ∂V ∂r • ∂ ∂r p • ∂ ∂r E(t) • ∂ ∂r V + O(ǫ 7 ). (16b)
If one needs to know the averaged Hamiltonian and/or the system of coordinates at the order ǫ N , one needs to truncate the O(ǫ N +1 ) terms of Eqs. ( 14) and/or (16), respectively. The expressions (16) are used below to reconstruct the trajectories of the electrons from the trajectories of the guiding centers.

In summary, the hierarchy of models is composed of a time-independent Hamiltonian H m [see Eqs. (15)] and a transformation from the electron coordinates to the guiding centers:

Φ n : (r, p) → (r, p), ( 17
)
whose truncated expressions are given by the truncations of Eqs. (16). The orders m and n refer to the order of the Hamiltonian and the transformation, respectively, after truncation of the perturbative expansion. We refer to the model as the n-th order guiding-center model, where n ≥ m, and we show the set of (m, n) in Table I.

G n = (H m , Φ n ), Φ 2 Φ 3 Φ 4 Φ 5 Φ 6 Φ 7 H 2 G 2 G 3 G 4 • • • H 5 • • • G 5 G 6 • H 7 • • • • • G 7

Links with the Kramers-Henneberger potential

At order ǫ 5 , our model is linked to the Kramers-Henneberger (KH) treatment of the motion of a charged particle in an external time-periodic electric field [43][44][45]. In a nutshell, the classical KH theory amounts to performing a canonical Lie transform generated by

S (2) = ǫ r • A(t) - ǫ 2 ω 2 p • E(t),
on Hamiltonian (8), where the first term is used for moving into the velocity gauge. The resulting Hamiltonian becomes exactly

H KH = k + ǫ |p| 2 2 + V r + ǫ 2 ω 2 E(t) + ǫ 3 E 2 (t) 2ω 2 .
The term of order ǫ 3 in the Hamiltonian H KH can be easily removed by performing an additional transformation. In KH theory, the remaining time-dependence in the potential is removed by an integral over time of the effective potential. For instance, if we expand the KH effective potential up to order ǫ 7 and we average over t ∈ [0, T ], it becomes

V KH = V + ǫ 4 U p ω 2 (ξ 2 + 1) ∂ 2 V ∂x 2 + ξ 2 ∂ 2 V ∂y 2 + O(ǫ 8 ).
In this framework, the effective potential always only depends on the position variables (and not on the momenta), and there is no contribution at order ǫ 6 , contrary to our derivations. The origin of this discrepancy is that performing an averaging using an integral over the fast timescales is only correct to the lowest order (here ǫ 4 ), but it fails at higher orders. One needs to perform canonical changes of coordinates to properly average the fast motions. Expressions beyond order ǫ 4 and results obtained using these higher orders are therefore incorrect.

We hereby take the opportunity to reiterate the advantage of using canonical Lie transforms in the reduction procedure: Since these transformations are invertible and their inverse can be algebraically computed, information on the original system [as described by Hamiltonian (1)] can be fully recovered using the dynamics of the reduced Hamiltonians. The model G 5 contains more information than what is provided by KH theory. In particular, the KH theory does not provide Φ 5 , and as a consequence, we are not able to reconstruct consistently the trajectory from H KH .

III. ANALYSIS OF THE GUIDING-CENTER MODELS

In this section, we analyze the different guiding-center models, composed of a time-independent Hamiltonian H m and a canonical transformation Φ n . In what follows, we restrict the analysis to a simplified potential for the atoms, namely the soft-Coulomb potential [46,47]:

V (r) = - 1 |r| 2 + 1 . (18)
A. Comparison between Hamiltonian (1) and the reduced models

Figure 1 shows typical electron trajectories (dark blue curve) of Hamiltonian (1), the guiding-center trajectories of G 2 (solid light blue curve) and G 5 (solid red curve) and the associated reconstructed trajectories (dashed curves), for d = 2, I = 10 14 W • cm -2 , ω = 0.05 a.u. Different ellipticities are considered: ξ = 0, ξ = 0.5 and ξ = 1. The guiding-center trajectories (r(t), p(t)) are computed by solving the forward and backward equations of motion of the corresponding guiding-center Hamiltonians (15) with initial conditions Φ n (r(t 0 ), p(t 0 )), where t 0 is chosen such that |r(t 0 )| > 50 a.u. (∼ 2E 0 /ω 2 ). We observe that the light blue and red solid curves guide the oscillating dark blue curves. Therefore, the electron oscillates around a guiding center. For instance, in Figs. 1(c)-(d), we see that the trajectory ionizes if the guiding-center motion is unbounded, and it returns to the ionic core if the guidingcenter returns to the core. Moreover, we observe a qualitative agreement between the electron trajectory and the reconstructed trajectory using the models G 2 and G 5 . In addition, we observe an overlap almost everywhere between the dark blue curve and the red dashed curve, a signature of a very good quantitative agreement between the electron trajectory and the reconstructed trajectory of the model G 5 . For the G 2 , the overlap is mainly observed far from the ionic core, i.e., for shorter integration times (less than 10T ). Below, we provide more thorough analyses to see if and when this agreement between the reduced models and the true trajectories holds. We consider the case d = 1 for clarity. In what follows, the parameters are ξ = 0, I = 10 14 W • cm -2 , and ω = 0.05.

As a consequence of d = 1, the electron and guidingcenter phase-space coordinates are reduced to (r, p) = (xx, px) and (r, p) = (xx, px), respectively. Looking at longer trajectories as it is done in Fig. 2, we observe mul- Lower panel: The dark blue curve is the guiding-center energy H5(Φ5(x(t), p(t))). The red curves are the guiding-center energy of G5, given by H5(Φ5(x(t0), p(t0))), respectively. The horizontal black line is E = 0. Here x and E are in atomic units.

tiple returns of the electron to the ionic core. The upper panel of Fig. 2 shows a typical trajectory (dark blue curve) of Hamiltonian (1), and the guiding-center trajectory for G 5 (red curve) for every interval of time when the electron is far from the ionic core. The guiding-center trajectory is solution of the forward and backward equations of motion of Hamiltonian (15b), with initial conditions Φ 5 (x(t 0 ), p(t 0 )), such that |x(t 0 )| > 50 a.u. In the lower panel of Fig. 2, the dark blue curve is the guidingcenter energy H 5 (Φ 5 (x(t), p(t))), i.e., at each time, the transformation Φ 5 is performed on the electron phasespace coordinates, and its associated energy H 5 is computed. The red curves are the guiding-center energy of Hamiltonian (15b) for initial conditions Φ 5 (x(t 0 ), p(t 0 )), i.e., the energy of the guiding center of G 5 , which is conserved. In Fig. 2, we observe that the guiding center reproduces well the mean trajectory of the electron for several time intervals when the electron is far from the ionic core, in a similar way as it was observed in Fig. 1. As a consequence, the guiding-center energy of the electron H 5 (Φ 5 (x(t), p(t))) is approximately conserved in a piece-wise manner in time. However, we notice that the energy strongly varies during close encounters between the electron and its ionic core. In addition, once the electron has undergone a close encounter, the guiding-center energy of the electron jumps to another energy level.

These observations on the reconstructed trajectories and on the guiding-center energy lead us to consider two In what follows, we use these two methods to test the validity and the benefits of the reduced models. Using method (i), the reconstructed trajectories of the model must be close to the true electron trajectory for the models to be relevant. This method is employed in Fig. 3. Using method (ii), by definition, the guiding-center energy of the electron H m (Φ n (r(t), p(t))) must be conserved up to some order for the reduced models to be relevant. This method is employed in Figs. 4 and5. In addition, we use these tests to compare the reduced models and to provide some guidelines on which models should be used for practical purposes.

Reconstructed trajectories

Figure 3 shows the most probable distance error δx n (t) between the reconstructed trajectories of the model G n and the electron trajectory, such that

δx n (t) = |Π(Φ -1 n (x(t), p(t))) -x(t)|, (19)
where x(t) is the trajectory of Hamiltonian (1), x(t) is the guiding-center trajectory of Hamiltonian H m with initial condition Φ n (x(0), p(0)), and Π is the projection from phase-space onto the position component, i.e., Π(x, p) = x. The most probable distance error is determined using the maximum of the kernel density estimation [48] of the distance error as a function of t/T . Specifically, it is determined in two steps for a fixed t/T : First, we compute the kernel density estimation of our data, then, we locate its maximum. The equations of motion for x and p for the models G 2 = (H 2 , Φ 2 ), G 3 = (H 2 , Φ 3 ) and G 4 = (H 2 , Φ 4 ) are the same, given that the Hamiltonian are the same. The differences between these models come from the change of coordinates, as taken into account in the determination of the initial conditions of the guiding-center trajectory and in the reconstruction of the trajectory from the guiding-center phase-space coordinates. This is also the case when comparing the models G 5 = (H 5 , Φ 5 ) and G 6 = (H 5 , Φ 6 ).

The distance error between the electron trajectories and the reconstructed trajectories is increasing for increasing time. Moreover, the distance error is increasing faster for n = {2, 3, 4} than for n = {5, 6, 7}. As a consequence, at t ≈ 100T , the distance errors δx n for n = {2, 3, 4} are two orders of magnitude greater than the ones for n = {5, 6, 7}. More quantitatively, the most probable distance error δx n scales as

δx n ∝ |t| αn .
We observe that α n ≈ 2.1 for n = {2, 3, 4}, and α n ≈ 1.1 for n = {5, 6, 7}.

Counter-intuitively, we observe no significant quantitative improvements, neither between the models G 2 , G 3 and G 4 , nor between the models G 5 , G 6 and G 7 . Hence, far from the ionic core, the corrective terms in the change of coordinates Φ 3 and Φ 4 , are negligible (at least for the chosen parameters), and the models G 2 , G 3 and G 4 provide similar results. In the same way, the corrective terms in Hamiltonian H 7 , compared with H 5 , are negligible, as well as the corrective terms in the change of coordinates Φ 6 and Φ 7 .

Guiding-center energy

We complement the analysis of the trajectories by looking at a specific property of the reduced models, namely the conservation of energy. In Fig. 4, an ensemble of trajectories is initiated such that the initial velocity of the guiding center G 5 is zero, and the initial position of the guiding center G 5 is normally distributed with a mean value 100 a.u. and a standard deviation 5 a.u. The initial laser phase is uniformly distributed φ ∈ [0, 2π]. The distribution in the lower panel of Fig. 4 represents the distribution of the guiding-center energy H 5 (Φ 5 (x(t), p(t))) as a function of t/T . In the upper panel, the dark blue curve is a typical electron trajectory in the ensemble. The red curve is the guiding-center energy H 5 (Φ 5 (x(t), p(t))) of the electron trajectory corresponding to the dark blue curve in the upper panel. For t < 7T , we observe in the lower panel that the distribution is peaked around the initial guiding-center energy of the electron. During this time, the electrons are far from the ionic core. At t ∼ 8T , the electrons get close to the ionic core, and the guiding-center energy distribution starts to spread out. When the electrons are close to the ionic core, their dynamics is highly nonlinear due to the competition between the strong laser and Coulomb fields, and the energy curve of a single trajectory in the lower panel starts varying significantly. In the meantime, the inset in the upper panel shows that the approximate trajectory of the reduced model no longer reproduces the electron trajectory. For t > 9T , the electron is far from the ionic core, and the red curve in the lower panel stops varying. It means that the model G 5 is again relevant, but for a different energy level than the initial energy. The arrow indicates the jump of the guiding-center energy after the close encounter with the ionic core. Close encounters with the ionic core are short time processes, and therefore cannot be averaged in time. It is expected that the fast-time average we perform fails to describe the various energy exchanges happening on these short time scales.

Figure 5 shows the most probable relative energy error δe n (t) for the models G n as a function of the distance between the electron and the ionic core |x(t)|, such that

δe n (t) = H m (Φ n (x(t), p(t))) -H m (Φ n (x(0), p(0))) H m (Φ n (x(0), p(0))) , (20)
where (x(t), p(t)) are the electron phase-space coordinates at time t. The most probable energy error is the maximum of the kernel density estimation [48] of the energy error. It is determined using the same technique as for computing the most probable distance error (see Fig. 3). The initial conditions are the same as in Fig. 3, and the integration is stopped when the electron reaches x = 1 a.u.

As expected, we observe that the energy error δe n increases when |x| decreases, i.e., as the electron approaches the ionic core. Far away from the ionic core, we observe that the most probable energy error δe n scales as δe n ∝ |x| -βn , with β n ≈ 3.0 for n = {2, 3, 4}, and β n ≈ 6.5 for n = {5, 6, 7}. As in Sec. III A 1, for we observe no significant quantitative improvements, neither between the models G 2 , G 3 and G 4 , nor between the models G 5 , G 6 and G 7 .

A cross-over is observed between all the models when the electron reaches ∼ 35 a.u. In particular, from Fig. 5, we observe that G 2 = (H 2 , Φ 2 ) gives the smallest energy errors among the reduced models close to the ionic core. The main reason is that H 2 and Φ 2 do no contain derivatives of the potential. As such, G 2 constitutes the most robust model among the hierarchy. Far from the ionic core (> 35 a.u.), the higher order models provide a better quantitative agreement with the electron trajectories as shown in Figs. 3 and5. The efficiency of the higher-order models appears far from the ionic core, around 35 a.u.

Discussion

The relative failure of the higher order models close to the ionic core can be understood by looking at the magnitude of the corrective terms in the Hamiltonians H m and in the changes of coordinates Φ n . Figure 6 shows h 1 = V (x), h 5 = (U p /ω 2 )V ′′ (x) and h 7 = (U p /ω 4 )V ′′ (x) 2 as a function of the distance between the guiding center and the ionic core. The term h n corresponds to the timeindependent term in Hamiltonian (14) of order ǫ n . We observe an overlap between these terms, where h 1 ∼ h 5 for |x| ≈ x c . Far from the ionic core, we approximate the soft-Coulomb potential by a hard-Coulomb potential V (x) ≈ -1/|x| and we compute h n explicitly. We deduce that h 1 ∼ h 5 for x ≈ x c with

x c ∼ E 0 /ω 2 ,
which is approximately equal to 21 a.u. for I = 10 14 W • cm -2 and ω = 0.05 a.u. Therefore, for |x| ≫ x c , the terms are ordered such that h 1 > h 5 > h 7 . In this case, the series in the perturbative expansion of Hamiltonian (14) and the change of coordinates (16) are likely converging, and the models G n for n > 2 are relevant for the guiding-center dynamics. For |x| ≪ x c , the terms are ordered such that h 1 < h 5 < h 7 . In this case, the series in the perturbative expansion of Hamiltonian (14) and the change of coordinates (16) are likely diverging, and the models G n for n > 2 are no longer relevant for the dynamics of Hamiltonian (1).

B. Guiding-center phase-space dynamics 1. For H2

The first (non-trivial) element of the hierarchy is G 2 = (H 2 , Φ 2 ). This model was identified above as the most robust one in the hierarchy for the analysis of the guidingcenter dynamics, since its guiding-center energy error is lower than for any other models close to the ionic core. We notice that this is the only reduced Hamiltonian which does not depend on the parameters of the laser field. The dependence on the laser field is in the change of variables Φ 2 . Moreover, if the potential is rotationally invariant, as is the case for atoms, the resulting Hamiltonian is integrable since the angular momentum is conserved in addition to the Hamiltonian.

The change of variables is exactly given by r = r -ǫ 2 E(t)/ω 2 , p = p -ǫA(t).

What is particularly convenient with this guiding-center model is that the potential is taken into account in the Hamiltonian and the electric field in the change of variables.

Ionization occurs if and only if the energy of the guiding center E = H 2 (r(t), p(t)) is positive. Otherwise the motion of the electron is bounded since the guiding center moves on a quasi-periodic orbit. The laser parameters have no influence on the motion of the guiding centers (and this holds up to the fourth-order model). They only influence how the electron swirls around the quasiperiodic orbit.

We consider the case when the guiding-center is far from the ionic core and as a consequence we can approximate V (r) ≈ -1/|r|. Thus, the typical guiding-center trajectory is on a Kepler orbit for E < 0, where E is the energy of the orbit, as it is the case for the guidingcenter trajectories in Figs. 1(a-c) and (f). One of the particularities of these orbits is that the radial momentum p • r/|r| vanishes twice in a revolution cycle: Once when the guiding-center trajectory is at the perihelion r (minimum distance from the ionic core) and d(p•r/|r|)/dt > 0, and once when the guiding-center trajectory is at the aphelion R (maximum distance from the ionic core) and d(p • r/|r|)/dt < 0. The aphelion and the perihelion are such that R + r = 1/|E|, imposing that for a given energy, the larger the aphelion, the smaller the perihelion, i.e., the closer the electron gets to the ionic core.

For H5 and H7

When going to higher-order models, the Hamiltonian H 2 gets perturbed by h 5 (and h 7 ). As a consequence, since the perturbation mainly affects the trajectories that pass close to the ionic core, i.e., r < x c ∼ E 0 /ω 2 , we shall see that the most perturbed trajectories in the higherorder models are the ones with a large aphelion 1/|E| -E 0 /ω 2 < R < 1/|E|.

The fact that the energy of the guiding centers is conserved is a property which is preserved by construction of the reduction procedure: E = H 5 (r(t), p(t)) or E = H 7 (r(t), p(t)) is conserved in time. Consequently, for d = 2, the dimension of phase space is reduced from 5 to 3 + 1. Here 3 + 1 means that phase space is foliated by constant energy surfaces of dimension 3. The advantage is that one can visualize the dynamics using Poincaré sections. Figure 7 shows the Poincaré sections r • p = 0 and d(p • r/|r|)/dt < 0 for E = -0.01 a.u. and the soft-Coulomb potential (18), where p •r/|r| = d|r|/dt is the radial momentum of the electron. This Poincaré section corresponds to the position of the guiding-center when it turns back towards the ionic core. The region close to the ionic core, i.e., around |r| ∼ E 0 /ω 2 (grey areas), is not relevant since the reduction procedure is only valid far away from the ionic core.

For linear polarization (ξ = 0), in Figs. 7(b) and (f), we observe two distinct dynamical behaviors of the guidingcenter trajectories. These figures display a chaotic layer far from the ionic core around 1/|E| -E 0 /ω 2 < |r| < 1/|E|. This ring corresponds to trajectories with a small perihelion and a large aphelion, that come close to the ionic core |r| < E 0 /ω 2 . These trajectories are the most affected by the perturbation h 5 in the Hamiltonian H 5 according to our discussion in Sec. III A 3, and typically correspond to the trajectories for which the electron comes back to the ionic core, as it is depicted in Figs. 1(a),(c) and (f). The width of this ring is of order E 0 /ω 2 . Secondly, we observe a regular region for 1/2|E| < |r| < 1/|E| -E 0 /ω 2 . This region corresponds to trajectories with a perihelion greater than E 0 /ω 2 , of the same order of their aphelion. These are the trajectories least affected by the perturbations h 5 and h 7 . As a consequence, we observe orbits that are mostly preserved from the unperturbed Hamiltonian H 2 , and typically correspond to the trajectories for which the electron stays far from the ionic core, as it is depicted in Fig. 1(b). Also, we observe two elliptic islands for (x, ȳ) ∼ (±1/2|E|, 0). These islands correspond to nearly circular guiding-center orbits, which become stable with the coupling with the electric field encapsulated in the effective potential of H 5 or H 7 . Finally, the Poincaré sections for H 5 and H 7 look similar. This is another indication that h 7 does not bring enough perturbation compared to h 5 for these parameters and far from the ionic core.

For elliptical polarization (ξ = 0.5), in Figs. 7(c) and (g), the observations are similar to the linear case, which reinforces the generality of the discussion above. However, for circular polarization (ξ = 1), in Figs. 7(d) and (h), we no longer observe chaotic behavior in the guiding-center dynamics. Indeed, Hamiltonians H 5 and H 7 are rotationally invariant, and as a consequence, the guiding-center angular momentum is conserved. There- fore, the dimension of phase space is reduced from 3 + 1 to 2 + 2 and the system is integrable. We observe that the elliptic islands we observed for ξ = 0 and 0.5 around (x, ȳ) ∼ (±1/2|E|, 0) are no longer present in the circular polarization case, as a consequence of the rotational invariance.

In the previous sections we have seen that for the parameters I = 10 14 W • cm -2 and ω = 0.05 a.u., the models H 5 and H 7 are almost equivalent since h 5 is several order of magnitude higher than h 7 in the region where the higher order models are relevant, |r| > E 0 /ω 2 (see Fig. 6). This is verified by comparing the Poincaré sections of Figs. 7(b-d) and Figs. 7(f-h). However, in Figs. 7 (e)-(i), we observe that when ω = 0.025 a.u. so that the quiver radius is of the same order as the distance between the guiding-center and the ionic core (E 0 /ω 2 ∼ 85 a.u.), the dynamics between Hamiltonians H 5 and H 7 differs significantly. The reduced models are not relevant in these regions. Hence, the reduced models G 5 and G 7 are significantly different when the characteristic distance between the guiding center and the ionic core is the same as the quiver radius (at least for this range of parameters). Similar observations would have been made for ω = 0.05 a.u. by lowering the guidingcenter energy E.

Conclusions

We have derived a hierarchy of reduced models G n for the guiding-centers dynamics of the electron interacting with the combined strong laser and Coulomb fields.

The reduced models G n are composed of an averaged Hamiltonian H m governing the guiding-center dynamics [Eqs. (15)] and a transformation Φ n which maps the electron phase-space coordinates onto the guiding-center phase-space coordinates [Eqs. (17)]. As a rule of thumb, these models are relevant when the electron is relatively far away from the ionic core (typically when its distance from the core exceeds one quiver radius), which happens in a piece-wise manner in time. The models do not describe the short events when the electron recollides with the ionic core.

We have singled out two models G 2 and G 5 : The first model provides the leading behavior of the trajectories and is the most tractable one due to its simplicity. In order to improve the quantitative agreement, a higherorder model such as G 5 has to be used.

All these models allow the distinction between direct ionizations and rescattering with the ionic core. This is very useful when the photoelectron momentum distributions are analyzed for imaging the target. In particular, we were able to define an energy of the electron far away from the core for this time-dependent system. The rescattering events can be seen as jumps in energy as a result of the transfer of energy from the ionic core to the electron.

FIG. 1 :

 1 FIG. 1: Typical electron trajectories for I = 10 14 W • cm -2 , d = 2 and ω = 0.05 a.u. in the polarization plane (x, y). The ellipticities are: (a) ξ = 0, (b-e) ξ = 0.5, (f) ξ = 1. The dark blue curve is the electron trajectory of Hamiltonian (1). The light blue and red curves are the trajectories of the models G2 and G5, respectively, with initial conditions far from the ionic core. For each model, the solid and dashed curves are the guiding-center and the approximate trajectories, respectively. The lightly shaded circle is the position of the ionic core at the origin, and the black circle surrounding the origin is |r| = 15 a.u. All quantities are in atomic units.

=

  z -{S, z} + 1 2 {S, {S, z}} + . . . , where L S is the Liouville operator defined by L S F = {S, F }. Canonical Lie transforms have several properties: (i) F (exp(L S )z) = exp(L S )F (z), which comes from the Leibniz rule, (ii) {exp(L S )F, exp(L S )G} = exp(L S ){F, G}, which comes from the Jacobi identity and the antisymmetry of the Poisson bracket,

FIG. 2 :

 2 FIG. 2: Typical electron trajectory of Hamiltonian (1) for I = 10 14 W • cm -2 , d = 1, ξ = 0 and ω = 0.05, and the guiding-center trajectory of the model G5, for multiple initial conditions Φ(x(t0), p(t0)) such that |x(t0)| > 50 a.u., with forward and backward integration of the equations of motion of H5. The grey areas are where the guiding-center position is |x| < 35 a.u. Upper panel: Dark blue and red curves are the electron and the guiding-center trajectory, respectively.Lower panel: The dark blue curve is the guiding-center energy H5(Φ5(x(t), p(t))). The red curves are the guiding-center energy of G5, given by H5(Φ5(x(t0), p(t0))), respectively. The horizontal black line is E = 0. Here x and E are in atomic units.

FIG. 3 :

 3 FIG. 3: Most probable distance error δxn [see Eq. (19)] as a function of time t per laser cycle T (in log-log scale) for I = 10 14 W • cm -2 , d = 1, ξ = 0 and ω = 0.05 a.u. The electrons are initialized such that the initial velocity of the guiding center of G2 is zero, and its initial position is normally distributed, with a mean value 1000 a.u. and a standard deviation 5 a.u. The initial laser phase is uniformly distributed φ ∈ [0, 2π]. The dashed lines are the linear fit curves. The background color is the mean value of the distance of the electron as a function of t/T . Here δxn is in atomic units.

FIG. 4 :

 4 FIG. 4: Electron trajectory for I = 10 14 W • cm -2 , d = 1 and ω = 0.05 a.u. The gray areas are where the guiding-center position is |x| < 35 a.u. Upper panel: The dark blue curve is the electron trajectory of Hamiltonian (1) as a function of time t per laser cycle T . The solid and dashed red curves are the guiding-center and the approximate trajectory for G5, respectively, with initial condition t0 = 0. The inset is a zoom of the region around the grey area. Lower panel: Logarithm of the distribution of the guiding-center energy H5(Φ5(x(t), p(t))) as a function of t/T . The guiding-centers are initialized with a normal distribution with mean value 100 a.u. and standard deviation 5a.u., with zero-velocity and a uniformly distributed initial laser phase φ ∈ [0, 2π]. The red curve is the guidingcenter energy H5(Φ5(x(t), p(t))) of the dark blue curve in the upper panel. Here x and E are in atomic units.

FIG. 5 :

 5 FIG. 5: Most probable energy error δen [see Eq. (20)] as a function of the distance between the electron and the ionic core |x| (in log-log scale), for I = 10 14 W • cm -2 , d = 1 and ω = 0.05 a.u. The initial conditions are the same as in Fig. 3. The dashed lines are the linear approximation for |x| ∈ [1, 35] a.u. The inset shows a zoom of the curves. Here |x| is in atomic units.

FIG. 6 :

 6 FIG. 6: Representation of hn as a function of the distance between the guiding center and the ionic core |x| (in log-log scale) for I = 10 14 W • cm -2 , d = 1 and ω = 0.05. The term hn corresponds to the term of order O(ǫ n ) in the Hamiltonian (14), which are: h1 = V (x), h5 = (Up/ω 2 )V ′′ (x), and h7 = (Up/ω 4 )V ′′ (x) 2 . The dashed lines are the asymptotic behaviors of these terms for soft-Coulomb potential (18), which are for |x| ≫ 1: |h1| ∝ |x| -1 , |h5| ∝ |x| -3 and |h7| ∝ |x| -6 . The vertical line is xc = E0/ω 2 . Here hn and |x| are in atomic units.

FIG. 7 :

 7 FIG. 7: Poincaré sections r• p = 0 and d(p•r/|r|)/dt < 0 in the polarization plane (x, ȳ), for potential (18), I = 1×10 14 W.cm -2 and E = -0.01 a.u. The grey areas are |r| < E0/ω 2 . The Hamiltonians are: (a) H2, (b-e) H5, and (f-i) H7. The ellipiticities are: (b,f) ξ = 0, (c,e,g,i) ξ = 0.5, (d,h) ξ = 1. The frequencies are ω = 0.05 except for (d,h) ω = 0.025 (larger quiver radius). All axes are the same as for (a) unless stated otherwise. All quantities are in atomic units.

TABLE I :

 I The guiding-center models Gn = (Hm, Φn), where n ≥ m.

Acknowledgments

The project leading to this research has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-