Carrier–envelope phase effect in the yield of sequential ionization by an intense few-cycle laser pulse

The relative yield of highly charged atomic ions produced by a short (4–6 fs at FWHM) intense (1014–5 × 1018 W cm−2) laser pulse was investigated by numerical solution of the rate equations. We predict oscillations of the ion yield as a function of the absolute phase. A distinctive property of this phase dependence is that it can only be observed when at least two ions have comparable yields. It is shown that with currently available laser systems the effect should be experimentally detectable for various rare gas atoms: Xe, Kr, Ar and Ne.


Introduction
Double and multiple ionization of atoms and molecules is one of the most theoretically challenging strong-field phenomena. This process has been intensively investigated over the last two decades, both experimentally and theoretically. As a result of these investigations it has been found that two different physical scenarios underlie this phenomenon: sequential and nonsequential double (multiple) ionization (for reviews see [1][2][3][4]). In the process of sequential ionization, electrons are ionized sequentially and independently of one another.
In contrast to this process, electron-electron correlation is a necessary condition for nonsequential double (or multiple) ionization. The actual mechanism underlying this correlation has been the subject of extensive investigation and debate for several years. At this point it has been established that nonsequential double ionization occurs due to inelastic recollision of the ionized electron with the parent ion [5,6]. This mechanism can be visualized by the well-known three-step model. According to this model, (i) one electron tunnels out of the atomic potential and is driven by the laser field away from the residual atom. When the field changes its direction, (ii) the electron is driven back, and on recollision with the ion core, (iii) it gives part of its kinetic energy to a second bound electron, which is then freed. Finally, both are driven by the laser field. Note that the details of this picture significantly depend on the specific atom (see, e.g., [1,2]). Obviously, the rate of the nonsequential process is not the product of single ionization rates. It was first proposed in [7,8] that the nonsequential mechanism contributes to double ionization.
In contrast to nonsequential ionization, which can be observed at relatively low laser intensities, sequential ionization dominates at high intensities. Sequential processes also dominate in the case of ionization by a circularly polarized laser field. For elliptical polarization the electron must depart with nonzero initial velocity perpendicular to the main axis of the polarization ellipse in order to return to its parent ion. However, this initial velocity decreases the tunneling rate. As a result, the probability of the nonsequential double ionization decreases quickly with increasing ellipticity and for a circularly polarized field, the nonsequential mechanism is essentially eliminated. At present, sequential ionization is attracting particular interest due to the constant increase in the intensity of modern lasers. Moreover, the rapid progress of the experimental technique (particularly, the advent of the cold-target recoil-ionmomentum spectroscopy COLTRIMS or REMI) has not only allowed for measurements of the total ion yield, but also the momentum distributions of the product ions [9][10][11][12][13][14][15].

3
Strange as it may seem, at present, there is a lack of theoretical studies of sequential double and multiple ionization. This particularly holds for the ion momentum distributions and for ionization by short laser pulses. At first glance, there is a considerable body of work on sequential processes, such as [16][17][18][19][20][21][22][23][24][25][26]. However, the vast majority is concerned with only the total ion yield. Moreover, emphasis has been put on the description of the well-known 'knee' in the yield, which appears due to the nonsequential channel. Momentum distributions of multiply charged ions generated by sequential ionization of atoms were studied in the works [24][25][26], two of which were experimental and theoretical investigations [24,25] of the momentum distributions of multiply charged He, Ne and Ar ions up to the charge of Z = 8 generated by a long (200 fs at full-width at half-maximum (FWHM)) laser pulse with the maximum intensity of 5-7 × 10 15 W cm −2 . In the third work [26], the momentum distributions of Xe ions up to Xe 24+ were calculated for both relatively long and short laser pulses.
To the best of our knowledge, our previous work [26] is the only theoretical study of sequential multiple ionization by a short laser pulse, when the carrier-envelope phase (CEP) effects may be essential (for a review of these effects, see [27] and references therein). In [26], we predicted a substantial shift in the maximum (centroid) of the ion-momentum distribution along the laser polarization as a function of the CEP. This effect should be experimentally detectable with currently available laser systems even for relatively long pulses, i.e. 25-30 fs.
At first glance, it is much easier to measure the variation of ion yield with the absolute phase than to measure the phase-dependent shift of the momentum distributions. However, CEP effects are often not apparent in measurements of the total yields, because they can be washed out by the integration over momentum. Nevertheless, the effect of the absolute phase on the total yield of highly charged ions has not been studied so far. In this paper, we investigate the CEP dependence of the ion yield by solving the system of rate equations for different charge states. We perform a comprehensive analysis of the absolute phase effect for various atomic targets (Ne, Ar, Kr and Xe) over a wide range of laser intensities and formulate the conditions under which it can be observed.

Theoretical approach
In our calculations, we use a short linearly polarized laser pulse with a sinus-square envelope: where ω is the carrier frequency, ϕ is the absolute phase and τ L = (2π/ω)n p is the duration of the pulse; here, in turn, n p is the number of cycles within the pulse. The pulse duration at FWHM can be estimated as one half of the τ L : τ FWHM ≈ τ L /2. Assuming that ion creation is a sequential process (i.e. there is no interaction between the electrons), we come to the following system of rate equations: dr 0 dt = −w 0 r 0 , Here r Z ≡ r Z (t, ϕ) is the probability for an ion to have charge Z at time t, N is the maximum charge state that can be achieved at a given laser intensity and w Z ≡ w Z (F(t, ϕ), I Z ) is the 4 ionization rate for an ion with the charge Z and ionization potential I Z . The initial conditions for the system of rate equations, equations (2), are r 0 (0, ϕ) = 1 and r Z (0, ϕ) = 0 for Z = 1, 2, . . . , N . As the probability of finding a particle with some charge is equal to unity, the system (2) is constrained to the condition r 0 (t, ϕ) + r 1 (t, ϕ) + · · · + r N (t, ϕ) = 1.
In order to solve system (2), we need to know the ionization probabilities. According to simple estimates, the ionization of the neutral Xe atoms, as well as that of the Xe 1+ and Xe 2+ ions, occurs in the barrier-suppression regime, and, strictly speaking, one needs to use ionization rates suitable for barrier-suppression ionization; see, e.g., [28,29]. The formula of [28,29] reduces to the usual Perelomov-Popov-Terent'ev (PPT) [30] or Ammosov-Delone-Krainov (ADK) [31] rates in the tunneling limit, when the laser field is relatively weak. The extrapolation of the tunneling formula into the barrier-suppression regime overestimates the ionization rate as compared to the barrier-suppression one. However, a substantial difference between these two formulae (say, at the level of 50%) occurs only when the laser field is at least 7-8 times greater than the corresponding barrier-suppression value, which is not true in our case. Numerical simulations show that the ionization of neutral atoms, as well as that of singly and doubly charged ions, occurs mostly at field strengths no more than 2 times stronger in comparison with the barrier-suppression fields (see figures 2(a) and (d)). Moreover, ions with Z > 2 are ionized well before the intensity reaches the corresponding barrier-suppression value. Thus, barrier-suppression ionization plays a minor role in the production of highly charged ions, both for long [24] and short [26] laser pulses, and the PPT or ADK formulae for the ionization probability, w Z (t), are appropriate. Thus, the tunneling rate for a level with ionization potential I p and angular and magnetic quantum numbers equal to l and m, respectively, is given by where F a = (2I p ) 3/2 is the atomic field, n * = Z / (2I p ) is the effective quantum number and C 2 kl is the asymptotic coefficient of the atomic wave function; see [32]. We neglect the Keldysh parameter γ = ω 2I p /F in the exponent of equation (4) and omit the pre-exponential factor √ 3π/F, because the latter arises from the averaging over the laser period. The probabilities r Z (t, ϕ) can be calculated in various ways. In order to solve the system of the rate equations (2) numerically, it is reasonable to use the Gear method [33] or other methods adapted to the solution of stiff differential equations instead of the standard Runge-Kutta method. The fact that the system (2) is stiff becomes evident if we compare the ionization rates w k−1 and w k entering the same equation for the probability r k (t, ϕ). Usually these rates differ by one or even more than one order of magnitude. Another technique is based on the simulation procedure, which was developed in [26] for the derivation of ion momentum distributions. For high laser intensities, when multiply charged ions are available and, as a consequence, the system (2) consists of a large number of equations, the latter method is faster than the former one. We used both approaches in order to check consistency. The results obtained by both methods coincide with each other.
Bearing in mind potential measurements with widely available Ti:sapphire chirped pulse amplification (CPA) systems, we perform calculations for the wavelength of 790 nm. When solving the rate equations (2), we restrict ourselves to a maximal intensity of 5 × 10 18 W cm −2 . At this intensity the energy of an ionized electron in the final continuum state is of the order of its rest energy and the formulae for the relativistic ionization should be used in equations (2) instead of the rate (4).

Results
Let R Z (ϕ) ≡ r Z (t = τ L , ϕ), Z = 1, 2, . . . , N , be the yield of an ion with a charge Z at the end of the laser pulse, i.e. at t → ∞. In the following, we consider only these final yields of different charge states. We investigate how the yields R Z depend on the CEP ϕ. This dependence is shown in figure 1 for several charge states at two laser intensities. It is seen that the relative yields are oscillating functions of the absolute phase. The oscillation amplitude is rather complicated, not a monotonic function of the laser intensity. The minima and maxima of the yield curves correspond to ϕ = 0, π/2, π, 3π/2 and 2π. Moreover, if the yield of the ion with a certain charge Z has a maximum, the yield of the ion Z − 1 has its minimum at the same phase, and vice versa. If only two different charge states can be detected after the end of the pulse, this feature can be easily explained by equation (3). In the general case, it follows from the analytic formulae; see section 3.
The oscillations of the ion yield at a given intensity can be characterized by the difference between the yields at ϕ = 0 and ϕ = π/2: see figure 1. The quantity d may be negative or positive depending on the intensity and ion charge. This parameter can be considered as an absolute characteristic of the effect. It is useful to introduce the relative parameter as well: .
The parameters d and α are two independent characteristics. This is evident from the following example: at the intensity of 3.0 × 10 15 W cm −2 (see figure 1(a) figure 2 over the range of intensities 10 14 -5.0 × 10 18 W cm −2 for two different pulse durations. As expected, the shorter the laser pulse, the more pronounced the absolute phase effect. We first discuss the yields, see figures 2(a) and (d). For any given ion there exists an intensity region where the yield of this ion differs from zero. The left boundary of this interval is referred to as threshold intensity (for more details, see, e.g., [17], where the corresponding field strength for a long pulse was estimated analytically), whereas to the best of our knowledge, there is no common term for the right boundary. Let us call it the depletion intensity. It is seen from figures 2(a) and (d) that only two or, at some intensities, three yields differ from zero after the end of the short laser pulse. This is also true for any given time instant within the pulse, see [26].
The striking feature of the CEP effect can be revealed if we compare panels (a) and (b) (or (d) and (e)) of figure 2: the extrema of the differences d correspond to such intensities at which the ion yields are close to each other. If, on the other hand, the yields of all available ions are substantially different at a given intensity, the corresponding absolute values of d are small. In the case when only one charge state exists after the interaction with the laser pulse, there is no noticeable absolute phase effect. Indeed, long steps in figures 2(a) and (d) for Z = 8 and Z = 18 in the intensity ranges of 2-4.5 × 10 16 and 6.5-9 × 10 17 W cm −2 , respectively, give rise to the zero difference d; see figures 2(b) and (e). Such a situation occurs every time a subsequent atomic shell is completely ionized (in our case these are 5p 6 5s 2 and 4d 10 4p 6 4s 2 ). Hence, the presence of at least two ions with comparable yield is a necessary condition for oscillations with the CEP to be detectable. This feature distinguishes the phase effect, which we predict here for the sequential multiple ionization, from the other CEP effects known from the ionization by short laser pulses.
Intensity variation has a dramatic effect on the phase dependence of the ion yield: The maxima and minima change their places. This is clearly seen from figure 3(a), where the phase dependence of Xe 5+ yield is shown for several different intensities. As a consequence, the difference d changes its sign with increasing intensity; see figures 2(b) and (e). The evolution of d with intensity for Xe 5+ ion is also shown in figure 3(b), which was obtained by performing a cut of figure 2(b). It is seen that near the threshold intensity, d is positive and close to zero. It increases up to a maximum with increasing intensity. After the maximum, the difference decreases, passes through zero and becomes negative. Then d decreases further, while increasing in absolute value, and has its minimum at a certain intensity. Finally, when the intensity approaches the depletion value, d is negative and tends to zero. It should be noted that, according to equation (3), at any intensity the sum of all the values of d is equal to zero; see figures 2(b) and (e).
It is seen from figures 2(c) and (f) that, in contrast to the absolute characteristic, the relative quantity α has its maxima close to the threshold and depletion intensities of each charge state. The reason is that the ion yield varies faster with intensity than the difference d. It is important to keep in mind that the parameter α can be big when the ion yield is close to zero. For this reason, in figures 2(c) and (f) we show this parameter only if the corresponding relative ion yield is more than 0.05. It is seen, however, that even with such a restriction the relative characteristic can reach a value of 0.8 for shorter pulse at high intensities.
The dependence of the CEP effect on the atomic target is illustrated by figure 4, where the absolute characteristic d is shown for Ne, Ar and Kr ions; see figure 3(e) for Xe. It is seen that the heavier the atom, the more pronounced is the phase dependence of the ion yield. This feature will be explained below; see section 3.

Volume averaging
It is evident from this analysis that the phase dependence of ion yield is very sensitive to small variations in the laser intensity. This brings up the question: does the effect survive after focal averaging, which is inevitable in the experiment? In order to clarify whether the CEP effect can be observed in the worst possible case of 3D configuration, we integrate the yields R Z over all space. We assume that the laser intensity is a focused beam of Gaussian spatial profile and evaluate the volume integral by changing the integration variable from volume to intensity, see [17,34,35]. This integrating over isointensity shells greatly simplifies the calculations. The results of the averaging are shown in figures 5(a)-(c). It it seen that the effect is quite detectable, especially near the threshold intensities. Moreover, averaging mostly affects the negative values of d, reducing them approximately by one half, and the range of differences becomes asymmetric. These results could be expected. Indeed, for any charge state Z there is a distinct boundary between two intensity ranges, one of which corresponds to the positive and the other one to the negative values of d; see figures 2(b) and (e). Focal averaging leads to summation of the contributions from different intensities lower than a given peak value I 0 . These contributions begin to compensate each other only if I 0 belongs to the second intensity range; otherwise no compensation occurs due to focal averaging and the CEP effect survives.

Analytic model
Although numerical simulations have given an insight into the phase dependence of ion yield, it would be highly desirable to have an analytic model of the effect to gain a better understanding of the underlying physics. Such a model can be developed on the basis of the approximate solution of the rate equations (2), derived in [17]. This approximate solution, in turn, is based on the following formal solution of equations (2): Two approximations were performed in [17] to simplify equations (7): firstly, the probabilities, r k , in the intergands of equation (7) were approximated by their decaying edges: exp(− t 0 w k−1 (t , ϕ) dt ). This is well justified, if the charge states exist sequentially in time and the temporal overlap of the different states is low, i.e. the population transfer from Z = N to Z = N + 1 occurs quickly and all the Z = N population is transferred to Z = N + 1 before there is a significant Z = N + 2 population; see [26]. Our numerical analysis of the ionization dynamics for a short laser pulse shows that the aforementioned conditions are fulfilled provided the intensity is less than ≈ 5-7 × 10 17 W cm −2 . For example, at an intensity of 4 × 10 14 W cm −2 , Xe + ionizes with 99% probability within 0.8 of the optical cycle. Secondly, the ionization rates w k (t, ϕ) are assumed to be narrow square functions of time. This is also justified, because in the case of a short pulse every ionization event occurs on one half or on two adjacent halves of the laser period, i.e. on a small fraction of the whole pulse duration τ L . Thus, the rates w k (t, ϕ) are really narrow functions of time and, as a consequence, they may be approximated by squares of different heights.
After these approximations, the probabilities of different charge states after the end of the pulse read as (see [17]) . . . , where is the total ionization probability for an ion with charge Z , i.e. the total area under the rate curve, which can be calculated numerically. It is seen that the dependence of the yield R k , k = 1, . . . , N − 1, on the absolute phase is governed by the same dependence of the total rates W k and W k−1 . This contrasts with the case of momentum distribution, where the strong dependence on the CEP originates from various ionization pathways, which are realized at different values of the phase due to different temporal evolutions of the field; see [26]. The CEP dependence of the total ionization rate with high accuracy can be approximated by the following ansatz: where W 0 k = W k (0) and W π/2 k = W k (π/2). Note that sine-like ansatz is a natural one for the absolute phase dependence of the single charged ions. Ion yields calculated from equations (8) and (9) are shown in figure 1 by dashed lines. It is seen that predictions of the analytic model are in good quantitative agreement with the numerical results. The only exception is Xe 17+ ion at the intensity of 5 × 10 17 W cm −2 ; see figure 1(b): equations (8) overestimate the ion yield by a factor of two. This disagreement is caused by the fact that at high laser intensities the charge states are not separated in time and the approximation of the ion yield by its decaying edge is no longer valid, see above. However, even in this case, the model qualitatively describes the shape of the phase dependence.
In the following, let us consider the yield of the ion with the maximal charge Z = N . After the substitution of equation (9) the last equation of the system (8) reads as For most intensities and ion charges the difference (W 0 N − W π/2 0 ) is small compared to unity. This allows us to simplify equation (10) by expanding the exponent in a Taylor series: The product (W 0 (11) is the oscillation amplitude of the ion yield with the absolute phase, i.e. one half of the difference d: d/2. We notice that the value of this amplitude is governed by the interplay of the probability W 0 N −1 and the difference . If only one charge state N exists after the end of the laser pulse, the ionization probability W 0 N −1 is much greater than unity, and the oscillation amplitude tends to zero due to the factor of exp(−W 0 N −1 ). On the other hand, both the probabilities W 0 N −1 and W π/2 N −1 and therefore their difference approach zero in the vicinity of the threshold intensity. As a result, the oscillation amplitude is again negligibly small and no phase dependence can be detected. Assuming that the difference of the probabilities (W 0 N −1 − W π/2 N −1 ) is some fraction of W 0 N −1 , one easily finds that the amplitude in equation (11) has its maximum at W 0 N −1 = 1 when R N ≈ 0.63 and therefore the yields of at least two ions substantially differ from zero.
The greater the ionization potential, the sharper the phase dependence of the total ionization rate, i.e. the more substantial the difference between the total rates W 0 N −1 and W π/2 N −1 . For this reason, the CEP effect in the ion yield is more pronounced for heavier atoms than for lighter ones. For example, let us consider the ionization of Ne and Xe at the same intensity of 2.5 × 10 17 W cm −2 , when Ne 8+ and Xe 13+ ions dominate in the corresponding yields. The ionization potentials of the previous charge states, i.e. Ne 7+ and Xe 12+ , are equal to 239.1 and 294.0 eV, respectively. For the ionization of Ne 7+ the difference of the rates is (W 0 7 − W π/2 7 ) ≈ 0.03, whereas for Xe 12+ this difference is (W 0 12 − W π/2 12 ) ≈ 0.16. As a result, for Xe 13+ one has d = 0.06, which is six times greater than the same characteristic d = 0.01 for Ne 7+ . Thus equation (11) explains all the basic features of the numerical results obtained in section 2.

Conclusion
In conclusion, on the basis of the rate equations we investigate the relative ion yields of sequential ionization generated by a short (4-6 fs at FWHM) intense (10 14 -5 × 10 18 W cm −2 ) laser pulse. It is shown that the ion yields are oscillating functions of the absolute phase of the pulse, provided the pulse length is not more than 6 fs at FWHM. In contrast to other wellknown CEP effects, the phase dependence of the yield can be observed only if, at the end of the laser pulse, the yields of at least two ions are comparable with each other. Despite the strong intensity dependence, the effect survives focal averaging and should be experimentally detectable. In addition to the numerical results, we present an analytic model which explains the basic features of the phase dependence.