High-order above-threshold ionization beyond the electric dipole approximation

Photoelectron momentum distributions from strong-field ionization are calculated by numerical solution of the one-electron time-dependent Schrödinger equation for a model atom including effects beyond the electric dipole approximation. We focus on the high-energy electrons from rescattering and analyze their momentum component along the field propagation direction. We show that the boundary of the calculated momentum distribution is deformed in accordance with the classical three-step model including the beyond-dipole Lorentz force. In addition, the momentum distribution exhibits an asymmetry in the signal strengths of electrons emitted in the forward/backward directions. Taken together, the two non-dipole effects give rise to a considerable average forward momentum component of the order of 0.1 a.u. for realistic laser parameters.


Introduction
The interaction of strong laser pulses with atoms is a rich research area with a wide range of phenomena that involve large numbers of absorbed photons. Above-threshold ionization (ATI) is one of the early manifestations of multiphoton physics [1]. The name ATI indicates that in this process more photons than necessary to overcome the ionization threshold are absorbed. This explains the typical appearance of strong-field photoelectron energy spectra consisting of peaks separated by the photon energy of the ionizing laser beam. Over the years it was established that lowest-order perturbation theory in the light-matter interaction is usually not adequate to describe ATI correctly. For example, the ATI peaks at low electron energy are often not the most probable ones [2,3] and the peak positions can be shifted away [4] from the energy values expected on the basis of the simple energy-conservation rule involving the ionization potential of the atom and the number of absorbed photons, i.e. the multiphoton generalization of Einstein's photoeffect law. Important non-perturbative treatments of strong-field ionization are the strong-field approximation [5][6][7] and the numerical solution of the time-dependent Schrödinger equation (TDSE). The simplest possible non-perturbative approach, however, is the two-step model, where laser-induced electron emission is a sequence of (i) ionization and (ii) acceleration of the electron as a classical particle in the laser field [8,9]. In step (ii), the atomic potential is usually neglected.
One may say that the most striking non-perturbative aspect of ATI emerged with the discovery of high-order above-threshold ionization (HATI), namely the presence of a plateau in the energy spectrum at high photoelectron energies [10,11]. Photoelectrons with much higher energies than expected for 'direct electrons' described by the two-step model were measured and explained in terms of rescattering: electrons that are driven back to the parent ion may scatter elastically and subsequently undergo a second stage of acceleration in the laser field. The return of the electron to the parent ion is known as the recollision step. This adds a third step which is the reason that we may speak of the three-step model of HATI. It is important to note that the acceleration after recollision is the reason for the increased electron energies observed in HATI, so it could be more appropriate to speak about a four-step model, which, however, would not be consistent with the commonly used terminology. Analyzing the classical trajectories and neglecting the Coulomb force predicts the cut-off ≈10 U p for the electron energy, where w = ( ) U E 4 p 0 2 2 is known as the ponderomotive potential for a linearly polarized laser field with amplitude E 0 and frequency ω. (Here and in what follows, we use atomic units unless otherwise stated.) This classical model is analogous to the three-step model of high-order harmonic generation, which had already been proposed earlier [12,13]. In highorder harmonic generation, the electron recombines with the parent ion at the time of recollision resulting in emission of a single photon. In that context, the three-step model predicts the famous photon-energy cut-off at 3.17 U p +I p with the ionization potential I p .
Strong-field ionization is usually discussed in the electric dipole approximation, where it is assumed that the incident electric field is spatially homogeneous over the field-atom interaction region and furthermore the magnetic field in the laser pulse is neglected. Mathematically, the electric dipole approximation means that the exact vector potential ( ) A r t , of the incident pulse is replaced by a purely time-dependent function ( ) A t (assuming that we work in the Coulomb gauge and set the scalar potential to zero). In this approximation, the laser field cannot transfer any momentum to the center of mass of the atomic system. In other words, there is no radiation pressure in the electric dipole approximation. This implies that in the case of single ionization of an atom that was initially at rest, the photoelectron momentum and the momentum of the residual ion must add up to zero. Beyond the dipole approximation, the momentum of the absorbed photons is transferred to the particles and therefore in laserinduced ionization, an effect on the momentum distributions of electron and ion is expected. From a fundamental point of view, the theoretical treatment is complicated because the dynamics of the electron-ion two-particle system does not decouple into center-of-mass and relative motion in the nondipole regime. Since, however, the mass of the ion is much larger than that of the electron, it is a good approximation to assume that the center of mass is at rest and coincides with the ion position and to use the electron-ion relative coordinates as laboratory-frame coordinates. In this case, the photoelectron momentum distribution is the same as the momentum distribution calculated for the relative motion. In the electric dipole approximation, when a spherically symmetric system such as an atom is considered, the Hamiltonian describing the quantum mechanical dynamics of the electron coordinate is symmetric under interchange of the forward and backward direction, where 'forward' denotes the laser propagation direction. Thus, for ionization out of a spherically symmetric initial state, the average momentum of the photoelectrons will have a zero component along the laser propagation axis. While this conclusion holds in very good approximation for the laser parameters used in practical table-top experiments, the transfer of photon momentum to the photoelectrons has been observed in two recent experiments. Smeenk et al have used circularly polarized fields and have measured a small positive average momentum component in forward direction [14]. The positive shift is reproduced by a non-dipole version of the strong-field approximation [15][16][17]. A detailed analysis [18][19][20] shows that tunneling ionization contributes a forward momentum ( ) I c 3 p to the photoelectron momentum with c being the speed of light, which explains the findings in [17]. The major part of the forward shift in circularly polarized light, however, is directly proportional to the kinetic energy E e gained due to the acceleration in the laser field and can be written as E e /c. Ludwig et al [21] and Maurer et al [22], on the other hand, have shown that a negative shift, i.e. pointing against the propagation direction of the light, is found for the low-energy electrons from ionization in a linearly polarized field. As an explanation, the momentum change during rescattering in the Coulomb field of the parent ion has been given. The backward shift is reproduced also by numerical solution of the TDSE beyond the dipole approximation [23,24] and by the non-dipole formulation of the quantum trajectory-based Coulomb-corrected strong-field approximation [25].
In principle, the theory of strong-field dynamics including phenomena beyond the electric dipole approximation has already a considerable history [26][27][28][29][30][31][32][33][34][35][36]. However, it seems that most of the earlier work focused on high-order harmonic generation or considered quite high laser intensities in order to reveal relativistic effects. The non-dipole effects for direct (non-rescattered) electrons as well as the low-energy rescattered electrons have been analyzed by a number of groups in the last few years, as explained above. The high-energy rescattering plateau in photoelectron distributions has been investigated in [31,32] for high intensities beyond those used in typical table-top experiments.
The purpose of our contribution is the analysis of HATI in the regime of moderate laser intensities, restricting ourselves to terms of the order 1/c in the Hamiltonian. We present results from the numerical solution of the non-dipole TDSE and interpret them with the help of the non-dipole three-step model. Our central findings will be a deformation of the momentum distribution shape, which follows from classical trajectories, and an additional forward/backward asymmetry which cannot be explained without a more detailed modeling of the scattering step in the three-step model.
We restrict ourselves to a two-dimensional quantum mechanical simulation in the plane spanned by the electric field and the light propagation direction, because we expect only negligible non-dipole effects on the momentum distribution in the third direction.

TDSE with leading-order non-dipole corrections
We consider the ionization process of a single-active-electron atomic model in two dimensions under the influence of a linearly polarized plane-wave laser pulse of three cycle duration, fully represented by the electric field with η=t − z/c and frequency ω. The incident pulse travels in z-direction with the speed of light c=1/α≈137. The electric field E points along the x-axis of the coordinate system and the corresponding magnetic field =B e E c z points along the y-axis. Using the Coulomb gauge, we can set the scalar electromagnetic potential to zero and define a vector potential that is related to the physical fields by = B Aand = - ¶ E A t . In practice, we obtain the vector potential by numerical integration of the given electric field as We use the soft-core potential which has a Coulombic −1/r behavior at large r and reproduces the ionization potential I p ≈0.9 a.u. of helium. The ground state of this potential serves as initial state of the quantum mechanical propagation. Here, we assume that the ion mass is much larger than the electron mass and hence, the position and momentum of the electron can be approximated by the relative electron-ion coordinate r and the relative electron-ion momentum p because the center of mass almost coincides with the ion position [17]. For the same reason, the reduced mass is set equal to the electron mass, m=1. In principle, the dynamics of the system is determined by the non-relativistic Hamiltonian including the atom-field interaction with full spatial dependence of the electromagnetic field. However, for the laser parameters under consideration here, an expansion of the vector potential in 1/c can be applied and only leading-order corrections in 1/c are considered [24]. Inserting this expansion into equation (3) and neglecting terms proportional to 1/c 2 results in the Hamiltonian This form of the Hamiltonian has a disadvantage: the nondipole term contains a direct coupling between the z-component of the position operator r and the momentum operator p. Hence, even in the case of a vanishing potential the momentum p is not a conserved quantity. To overcome this drawback, we apply a unitary transformation and using its form invariance, we can read off the transformed HamiltonianH , (5), the transformation would be a simple gauge transformation, which is, in the limit 1/c→0, reminiscent of the transformation used in [38] to remove terms that do not explicitly contain the momentum operator. In our case, since position and momentum operators do not commute, the transformation results in a modification of the potential term in addition to the changes in the kinetic and interaction term. Using Hadamard's lemma leads to the following expression for the transformed Hamiltonian: The first term is familiar from the generalized Volkov states including leading-order non-dipole corrections derived in [20,27,29]. In the second term, the potential exhibits a position-and time-dependent shift z t c, which can be understood as a time-dependent shear of the potential.
The TDSE with the HamiltonianH is solved numerically using the standard Fourier split-operator method [39] with a time step of Δt=0.005 a.u. The numerical grid size is 209 a.u. in x-and 819 a.u. in z-direction with spacings of Δx=0.1 a.u. and Δz=0.2 a.u. The special form of the transformed Hamiltonian suggests a division of space into an inner interaction region and an outer asymptotic region so that the two regions are represented on two separate grids. In the inner part the complete HamiltonianH is implemented. In contrast, in the outer part the asymptotic tail of the potential is neglected and the electron is treated as a free particle interacting with the external electromagnetic field [40]. As the first term of the Hamiltonian is diagonal in momentum space, the time propagation in the outer region becomes trivial. The decomposition of the wave function into an inner and outer part is done by applying a mask function of 50 a.u. width on each side of the inner grid, optimized for electrons around » | | p 2.5 a.u., in every time step. What is absorbed from the inner region is transformed to momentum space and added to the wave function in the outer region. To speed up the time propagation and avoid the repeated computation of exponentials of the potential, we expand the shifted potential to first order in 1/c, i.e., The time-independent part of equation (8) is fully incorporated via exponential short-time propagators, while the small correction due to the second term is implemented as a Crank-Nicolson propagator. Taking into account higher-order terms of this expansion does not change the results presented below. After the end of the pulse the simulation is run for seven additional cycles to collect slow photoelectrons. The photoelectron momentum distribution ( ) p w is obtained at the final time from the asymptotic grid with a resolution of Δp x =0.015 a.u. and Δp z =0.004 a.u. Figure 1 shows the photoelectron momentum distribution for a pulse with 814 nm central wavelength corresponding to ω=0.056 a.u. and an intensity of I=8×10 14 W cm −2 . The momentum distribution has a typical shape consisting of a strong contribution of low-energy electrons concentrated along the p x -axis and a weaker contribution of high-energy electrons, which extend over a large range of momenta p z . In the energy region below the classical 2U p cut-off for nonscattered electrons (here at » | | p 2.6 a.u.

Results and discussion
x ) holographic interference stripes [41][42][43] overlaid with intracycle interferences are visible. The asymmetry between positive and negative p x is a familiar feature for very short few-cycle pulses [44]. Due to the shortness of the pulse only one dominating period of ionization times is selected and ATI peaks (i.e. peaks separated by the photon energy) are mostly avoided. The carrierenvelope phase of the pulse in equation (1) is such that the rescattering electrons with the highest energy are emitted with positive momenta p x and that their cut-off energy is close to that of a cw field. In the region of high-energy rescattered electrons, we clearly observe interference rings which arise from the superposition of signals from short and long rescattering trajectories [45,46]. The symmetry breaking in propagation direction (z-direction) is a non-dipole effect. It is too small to be directly visible in the full 2D spectrum. Therefore we show in figure 2 1D slices trough the distribution at different momenta p x . For the slices at p x =1, 2 and 3 a.u., we have averaged over an interval of Δp x =0.1 a.u. to suppress the strong oscillations from intracycle interference. In all energy regions we observe the following two major deviations from the dipole approximation: (i) the minima and maxima resulting from interferences are shifted in z-direction, i.e. their positions are no longer symmetric about the polarization axis; (ii) also the strengths of the peaks are different in forward and backward directions. In the following we refer to the second effect as asymmetry.
For small momenta (see figure 2(a)) the counter-intuitive shift of the central peak against the propagation direction discussed in [21,22,47] is found, in agreement with [23] and also the peak heights for emission in forward and backward directions are different [23]. Going to the higher momentum p x =2 a.u., this asymmetry becomes smaller and the shift of the center fringe becomes positive, in agreement with the simple classical estimation ( ) p c 2 x 2 [14,23]. In the intermediate energy region above the classical limit, but below the region of HATI, the maxima of the distribution are shifted in forward direction, but the asymmetry is inverted: in the slice at p x =3 a.u., the maximal signal in backward direction is higher than in forward direction, see figure 2(c).
We define a partial average that quantifies the momentum shift in propagation direction at fixed p x as ò ò This observable nearly follows the classical prediction ( ) p c 2 x 2 in the low and intermediate energy region, see figure 1(b). The oscillations below the classical 2U p cut-off are caused by the intra-and intercycle interferences. Deviations from the classical estimation can partially be attributed to the neglected tunneling process [18,19] and to the influence of the long-range coulomb potential [21,22].
In the high-energy region the 8-like structure with a cutoff close to the classical value of E≈10U p for positive momenta is clearly present in the momentum distribution of figure 1(a) [10,11]. The nearly circular ring structure is caused by interference between short and long rescattering trajectories [45,46]. In contrast to lower energies, a plateaulike structure appears in the partial average á ñ p z p x , see the range between p x =4 a.u. and p x =6 a.u. in figure 1(b). It differs strongly from the estimate ( ) p c 2 x 2 . Figure 2(d) shows that in this region the HATI signal exhibits a shift as well as an asymmetry with larger peak height in forward direction.
For a quantitative interpretation of the shifts in HATI, we extend the three-step model [11,13] beyond the dipole approximation. The three-step model consists of (i) laserinduced ionization, (ii) acceleration of the electron away and back to the parent ion and (iii) scattering from the parent ion with subsequent further acceleration in the laser field. The ionization step launches an electron with a velocity close to zero [19,31,32]. After the ionization takes place at time t i , the potential-free acceleration in the continuum can be described by Newton's equation, which gives rise to an electron For v/c=1 the equations can be solved perturbatively. If the initial velocity v 0 is chosen zero (which is the usual choice in the three-step model in the dipole approximation), we find that in leading order of 1/c only the motion in light propagation direction is modified by the non-dipole part of the Lorentz force. This drift motion causes the electron to miss the core and thus suppresses rescattering. For an efficient recollision process the electron has to start with an additional initial velocity = -| | v e v z 0 0 against the propagation direction of the light [27]. One finds that the required initial velocity is, in leading order, proportional to 1/c. For the resulting trajectory, one finds that the motion along the electric field axis is still not modified in order 1/c; only the motion along the light propagation is changed by terms of order 1/c. The associated ionization probability is slightly reduced compared to the dipole limit [27] due to to the non-zero initial velocity. But in contrast to [31,32], under the moderate laser parameters used throughout our work, the effect is only hardly seen. As we can always tune the initial velocity in propagation direction in a way that the electron returns to z=0, we only have to take the condition for exact return in polarization direction, x(t i )=x(t r ), into account to obtain the mapping between the ionization time t i and return time t r as in the dipole limit.
Using a classical Hamiltonian given by the expression as in equation (7), but without the potential term, we find that the potential-free trajectory starting from the position of the atom with initial velocity corresponds in leading order of 1/c to a conserved canonical momentum of = k - 0 . During rescattering, the electron feels the presence of the potential so that the dynamics is governed by the full Hamiltonian of equation (7) and the canonical momentum is changed from k to p. The latter is the final measurable momentum once we assume that no further interaction with the potential occurs. Under the assumption that the duration of the scattering process is short compared to period of the laser field, we can consider ( ) A t as constant during rescattering. Furthermore, since the rescattering takes place at z=0, the potential term - can be considered temporally constant. Therefore, the Hamiltonian is conserved during rescattering so that the momentum-dependent part before and after rescattering has the same value. For rescattering at time t r , we thus have In leading order in 1/c, the recollision energy on the righthand side of equation (10) is the same as in the dipole limit, depending only on the ionization and recollision time.
We define the magnitude of the corresponding velocity as = - r i r i . For fixed times, equation (10) selects a one-dimensional submanifold of momenta out of the complete two-dimensional space of possible momenta p. From a physical point of view, immediately after scattering the velocity of the electron lies on a circle of radius K(t r , t i ). Subsequently the electron is accelerated in the field for a second time and mapped to the final momenta. This acceleration process deforms the circle into an ellipse and, in leading order of 1/c, shifts its center by d = - . The principal axes are rotated by 45°relative to the xz-frame defined by the light polarization and propagation. If rescattering takes place at a time t r with < ( ) A t 0 x r , the center of the ellipse will be located in the region of positive p x ; in this case the semi-minor and semi-major axes are , respectively. By calculating all possible times and ellipses, we can determine the classical boundary of the momentum distribution. A schematic illustration of the resulting shape of the momentum distribution is given in figure 3. The deviation of the rescattering ellipse from a circle is quantified by the numerical eccentricity for an electric field with amplitude E 0 and frequency ω. Hence, the deformation of the momentum distribution shape due to non-dipole effects is predicted to increase with laser intensity and wavelength. In order to compare the TDSE results and the classical model, we regard the positions of the outermost interference maxima as the boundary of the TDSE momentum distribution. A numerically reliable value that quantifies the shifts of the TDSE interference peaks in the light propagation direction is given by The shift of the outermost ring is in perfect agreement with the classical prediction. Hence, the effect can be completely attributed to the acceleration of the electron in the electromagnetic field after rescattering. The gained additional forward momentum is equal to the energy gain after rescattering divided by the speed of light, which is evident from the last expression in equation (13) where p 2 2 is the final energy and K 2 /2 is the energy immediately after rescattering. This conclusion is plausible since both for classical light fields as well as photons, it is well known that every portion of light energy E traveling at light speed is accompanied by a momentum E/c. However, it is important to note that light energy absorbed before the rescattering event has no impact on the non-dipole shift of the boundary of the photoelectron distribution according to the three-step model.
Within each interference stripe, we find that the shift is a monotonic function of p x in the range p x >0. For small momenta p x we observe negative shifts that are bounded by at p x =0. In contrast to the outer ring, the classical model does not offer a simple procedure to interpret the inner rings, as for a given final momentum p long and short trajectories correspond to different rescattering velocities and therefore show different shifts. These trajectorydependent shifts are not directly visible in the distribution. Instead, the inner rings are caused by interference and strongly depend on the phase difference between the trajectories. Hence, to investigate these slightly smaller shifts, more advanced models such as the strong-field approximation in first Born-approximation have to be used [20,31,32].
As pointed out above, the photoelectron distribution exhibits not only a shifted and deformed shape, but also a forward/backward asymmetry in the peak heights. We quantify the asymmetry at each p x by taking the ratio of the signals of one interference maximum in forward direction and one interference maximum in backward direction, where each signal is obtained by integrating over one peak in p z -direction. The result is shown in figure 4(b). The asymmetry ratio is larger than one for the outer rings. It decreases while going to the inner rings and can reach values smaller than unity close to the center of the distribution. The ratio of the outer ring has a maximum of ≈1.158 located at p x ≈1.13 a.u. The combined effect of the shift and the asymmetry explains the plateau-like behavior in the partial average á ñ p z p x depicted in figure 1(b). Considering only the shift without asymmetry would yield a smaller partial average that matches with the full calculation only near the cut-off region, where the asymmetry has to vanish. Moving from the cut-off to smaller momenta p x , the decreasing shift is nearly compensated by the increasing asymmetry.
The shape of the plateau depends on the laser parameters as well as the atomic target. The first dependence is investigated in figure 5, where in panel (a) the intensity is scaled for a fixed wavelength of 814 nm and in panel (b) the wavelength is scaled for a fixed intensity of 8×10 14 W cm −2 . In all calculations we find the plateau-like structure. The height of the plateau increases with increasing intensity and wavelength. Additionally, the slope  over the plateau changes from negative values at small intensities or wavelengths to positive values at high intensities and wavelengths. We also find that the shifts extracted from TDSE are in perfect agreement with the classical model over the whole parameter range and that the dependence on the ring number is only weak as in the example depicted in figure 4(a). In contrast, the asymmetry drops faster from the boundary of the distribution towards small | | p z . Hence, for distributions that are wide in p z -direction as in the case of high wavelength and intensity, the influence of the shift dominates over that of the asymmetry, resulting in a positive slope over the plateau.

Conclusion
We have presented a numerical quantum mechanical investigation of HATI beyond the electric dipole approximation and an analytical analysis derived from the three-step model. To obtain photoelectron momentum distributions from the TDSE in a numerically convenient manner, we have used a unitary transformation that casts the Hamiltonian into a sum of purely momentum-dependent terms and purely position-dependent terms To our knowledge, such a form of the non-dipole Hamiltonian has not been used in any previous simulation. Our central result is the p x -dependent average momentum component in light propagation direction. In agreement with earlier studies, we find backward shifts at small energy and forward shifts at high energy. Additionally, we find that the p x dependence of the forward momentum in the high-energy rescattering regime exhibits a plateau-like structure at a considerable momentum of the order of 0.1 a.u. for the realistic laser parameters that we have chosen in the calculations. While the shift of the boundary of the momentum distribution can be explained by classical trajectories, there is an additional asymmetry in the forward/backward signal strengths. It appears that the explanation of this asymmetry must be sought in a non-isotropic rescattering event, which we plan to model by using appropriate elastic scattering cross sections in future work. Although all results shown in this paper have been calculated for short pulses, we have confirmed that the trajectory-based model yields very similar results for a cw monochromatic field. The non-dipole deformation of the momentum distribution shape becomes more pronounced with increasing laser intensity or wavelength. This has a clear interpretation in the classical model, which predicts that the rescattering region of the two-dimensional momentum distribution is composed of ellipses with eccentricity proportional to w ( ) E c 0 .