Frustrated tunneling dynamics in ultrashort laser pulses

In this subsection we review the SFA-FT theories presented in references [12, 13] in the language of Dyson equations. Consider a Hamiltonian H(t) that can be separated in two parts, H₀(t) and V(t), where the exact form of the time evolution operator associated to H₀ is known, U₀(t, t_i ). Then one can write down a recursive solution for the time evolution operator associated to the full Hamiltonian as:

$$\begin{equation}U(t,{t}_{i})={U}_{0}(t,{t}_{i})-i{\int }_{{t}_{i}}^{t}\mathrm{d}{t}^{\prime }\enspace U(t,{t}^{\prime })V({t}^{\prime }){U}_{0}({t}^{\prime },{t}_{i}).\end{equation} \tag{ 4 }$$

Using this Dyson equation for U, one can rewrite the exact amplitude for strong-field excitation, equation (1), as:

$$\begin{align}{a}_{n\ell m}& =\underset{t\to \infty }{\mathrm{lim}}-i{\int }_{{t}_{i}}^{t}\mathrm{d}{t}^{\prime }\\ & \quad {\times}\left\langle {\psi }_{n\ell m}(t)\right\vert U(t,{t}^{\prime }){V}_{L}({t}^{\prime }){U}_{A}({t}^{\prime },{t}_{i})\left\vert {\psi }_{0}({t}_{i})\right\rangle ,\end{align} \tag{ 5 }$$

with H₀ = H_A being the atomic Hamiltonian, and V(t) = V_L (t) being the light–matter interaction. Alternatively, the exact amplitude, equation (5), can be rewritten by an additional application of the Dyson equation, with the atomic potential, V(t) = V_A , as a second interaction to yield,

$$\begin{equation}\begin{aligned}{a}_{n\ell m}& =\underset{t\to \infty }{\mathrm{lim}}\left[-i{\int }_{{t}_{i}}^{t}\mathrm{d}{t}^{\prime }\left\langle {\psi }_{n\ell m}(t)\right\vert {U}_{L}(t,{t}^{\prime })\right.\\ & \quad {\times}{V}_{L}({t}^{\prime }){U}_{A}({t}^{\prime },{t}_{i})\left\vert {\psi }_{0}({t}_{i})\right\rangle \\ & \quad +{(-i)}^{2}{\int }_{{t}_{i}}^{t}\mathrm{d}{t}^{{\prime\prime}}{\int }_{{t}_{i}}^{{t}^{{\prime\prime}}}\mathrm{d}{t}^{\prime }\left\langle {\psi }_{n\ell m}(t)\right\vert U(t,{t}^{{\prime\prime}}){V}_{A}\\ & \quad {\times}{U}_{L}({t}^{{\prime\prime}},{t}^{\prime }){V}_{L}({t}^{\prime }){U}_{A}({t}^{\prime },{t}_{i})\left\vert {\psi }_{0}({t}_{i})\right\rangle \left.\right].\end{aligned}\end{equation} \tag{ 6 }$$

It should be noted that even if the recursive Dyson procedure is continued in this way, U will always appear in the last term of the expansion. Assuming that the exact form of U is not known, this implies that some approximation of the exact propagator must always be made in analytical calculations. In reference [13], Hu et al use equation (6) and approximate U by a propagator for dressed Rydberg states, U_nℓm . Formally, the FT process is then a second-order process that requires two interactions: (i) an excitation by the laser field at t' and (ii) a capture by the atomic potential at t'' > t'. This is because the first-order term with U_L vanishes in equation (6), due to quantum diffusion in the limit t → ∞ [13]. In this picture, the capture can happen at any time after the laser excitation, t'' > t', while in reality the population to Rydberg states via FT has to occur at a finite time before the end of the laser pulse. Less formally, first-order theories for FT can be developed by applying different approximations to U in equation (5) for different times based on physical intuition. For instance, one could approximate U by U_L when the field is strong, and then switch back to U_A when the field is weak. Provided that the laser field ends abruptly at finite time, t = t_f , Popruzhenko has proposed that the SFA momenta can be directly matched to Rydberg states using the Laplace–Runge–Lenz (LRL) vector [12]. In general, the 'end' of the laser pulse is not always possible to define uniquely, but for the ultrashort laser pulses of cos²-type that are considered in this work, defined in equation (3), the field ends at t_f = τ/2. This makes it possible to test the Popruzhenko SFA-FT theory with physical pulses and compare to Rydberg populations from exact TDSE simulations in hydrogen.

In order to relate the exact Dyson solution with the SFA theory at t_f , the exact amplitude for strong-field excitation in equation (1) can be written in terms of plane-wave momentum states using the identity, $\mathbb{1}=\int {\mathrm{d}}^{3}p\left\vert \mathbf{p}\right\rangle \left\langle \mathbf{p}\right\vert$, as

$$\begin{equation}\begin{aligned}{a}_{n\ell m}({t}_{f})=\int {\mathrm{d}}^{3}p\left\langle {\psi }_{n\ell m}({t}_{f})\vert \mathbf{p}\right\rangle a(\mathbf{p},{t}_{f}),\end{aligned}\end{equation} \tag{ 7 }$$

where

$$\begin{equation}a(\mathbf{p},{t}_{f})=-i{\int }_{{t}_{i}}^{{t}_{f}}\mathrm{d}{t}^{\prime }\left\langle \mathbf{p}\right\vert U({t}_{f},{t}^{\prime }){V}_{L}({t}^{\prime }){U}_{A}({t}^{\prime },{t}_{i})\left\vert {\psi }_{0}({t}_{i})\right\rangle .\end{equation} \tag{ 8 }$$

While equation (7) is exact and will approach a constant value for a_nℓm (t_f → ∞), it is worth to note that this formulation is not 'well behaved' due to $\left\vert \mathbf{p}\right\rangle$ not being an eigenbasis of H_A . As a result, the intermediate projection amplitude, a(p, t_f ), does not approach a constant value as t_f → ∞.

In the following, we have chosen to study the SFA-FT theory of reference [12] due its two main advantages: (i) the model is based on a first-order interaction with the laser field and (ii) all momentum variables are eliminated by substitution using the LRL vector.

From equations (7) and (8) it is evident that integration over both plane-wave momenta, p, and interaction times, t', must be performed in order to establish the SFA-FT amplitude. In saddle-point approaches, the order of integration is not important [16], and we first consider the momentum integration separately. Integration of momentum was omitted in prior work [12]. To proceed, we replace a(p, t_f ) in equation (7) by the SFA amplitude for photoionization

$$\begin{align}{a}_{\text{SFA}}(\mathbf{p},{t}_{f})& =-i{\int }_{{t}_{i}}^{{t}_{f}}\mathrm{d}{t}^{\prime }\\ & \quad {\times}\left\langle \mathbf{p}\right\vert {U}_{L}({t}_{f},{t}^{\prime }){V}_{L}({t}^{\prime }){U}_{A}({t}^{\prime },{t}_{i})\left\vert {\psi }_{0}({t}_{i})\right\rangle \end{align} \tag{ 9 }$$

and perform a saddle point integration over the momentum space,

$$\begin{equation}\begin{aligned}{a}_{n\ell m}({t}_{f})\approx \sum\limits _{{\mathbf{p}}_{s}}\left\langle {\psi }_{n\ell m}({t}_{f})\vert {\mathbf{p}}_{s}\right\rangle {a}_{\text{SFA}}({\mathbf{p}}_{s},{t}_{f}),\end{aligned}\end{equation} \tag{ 10 }$$

where the stationary momentum contributions are

$$\begin{equation}\begin{aligned}{a}_{\text{SFA}}({\mathbf{p}}_{s},{t}_{f})& \approx -i{\int }_{{t}_{i}}^{{t}_{f}}\mathrm{d}{t}^{\prime }{\left(\frac{2\pi i}{{t}_{f}-{t}^{\prime }}\right)}^{3/2}\\ & \quad {\times}\left\langle {\mathbf{p}}_{s}\right\vert {U}_{L}({t}_{f},{t}^{\prime }){V}_{L}({t}^{\prime }){U}_{A}({t}^{\prime },{t}_{i})\left\vert {\psi }_{0}({t}_{i})\right\rangle .\end{aligned}\end{equation} \tag{ 11 }$$

In writing equations (10) and (11), we have used the SFA action,

$$\begin{equation}{S}_{0}({t}^{\prime },\mathbf{p})={\int }_{{t}^{\prime }}^{{t}_{f}}\mathrm{d}{t}^{{\prime\prime}}\frac{1}{2}{[\mathbf{p}+\mathbf{A}({t}^{{\prime\prime}})]}^{2}-{I}_{p}({t}^{\prime }-{t}_{i}),\end{equation} \tag{ 12 }$$

and we have assumed that the momentum dependence on p of the Rydberg state, $\left\langle {\phi }_{n\ell m}\vert \mathbf{p}\right\rangle$, is small as compared to the SFA amplitude, a_SFA(p, t_f ). Rough estimations can be made using the width of the Rydberg state: Δp^Ryd ≈ 2/Δr^Ryd ≈ n⁻², while the SFA width is: ${\Delta}{p}^{\text{SFA}}\approx {({t}_{f}-{t}^{\prime })}^{-1/2}$. Physically, the window of capture by the Rydberg state in space should be small compared to the size of the returning SFA wave packet. The momentum projection, $\left\langle {\phi }_{n\ell m}({t}_{f})\vert {\mathbf{p}}_{s}\right\rangle$, can be approximated by its angular factor, $\left\langle {\phi }_{n\ell m}({t}_{f})\vert {\mathbf{p}}_{s}\right\rangle \approx {Y}_{\ell m}^{{\ast}}(\mathbf{p}/\sqrt{{p}^{2}})$ or evaluated using its analytical form, which is known for hydrogen [17]. A direct comparison of these prefactors will be shown in section 3.2. Clearly, equation (10) can be interpreted as FT because the electron goes from the ground state to a Rydberg state via the continuum states. The laser-driven dynamics is associated with a quantum diffusion factor, ${({t}_{f}-{t}^{\prime })}^{-3/2}$, that strongly damps long trajectories in the continuum. In SFA the condition of stationary momentum, p = p_s (t_f , t'), is that the electron returns to the atom,

$$\begin{equation}{\mathbf{r}}_{s}({t}_{f})={\int }_{{t}^{\prime }}^{{t}_{f}}\mathrm{d}\tau [{\mathbf{p}}_{s}+\mathbf{A}(\tau )]=({t}_{f}-{t}^{\prime }){\mathbf{p}}_{s}+{\Delta}G({t}_{f},{t}^{\prime })=0,\end{equation} \tag{ 13 }$$

which is a reasonable condition for a recapture process by the atom. Improvements to the recapture process on complex targets can be studied by adding a constant term, Δr_nℓ , on the right-hand side of equation (13). This implies a shift of the stationary momentum,

$$\begin{equation}{\mathbf{p}}_{s}=-{\Delta}G({t}_{f},{t}^{\prime })\to {\mathbf{p}}_{n\ell ,s}=(-{\Delta}G({t}_{f},{t}^{\prime })+\mathbf{\Delta }{\mathbf{r}}_{n\ell })/({t}_{f}-{t}^{\prime }),\end{equation} \tag{ 14 }$$

such that the final position, r_nℓ,s (t_f ) = Δr_nℓ , is reached within any given time window, t_f − t'. This ad hoc modification does not affect the quantum diffusion factor in equation (11), and it can be motivated by exterior centres of Rydberg states (for ℓ > 0) [13]. In the next subsection, we will show that the momentum substitution by the LRL-vector, proposed by Popruzhenko [12], leads to quasi-classical trajectories that do not return exactly to the atomic core at the end of the laser pulse.

In SFA the integral over interaction times, t' in equation (9), is evaluated using the saddle point method for which the action at time t' = t_s (p) is stationary by the tunneling condition,

$$\begin{equation}{\left[\mathbf{p}+\mathbf{A}({t}_{s})\right]}^{2}=-2{I}_{p}.\end{equation} \tag{ 15 }$$

The saddle point contributions are then weighted by the factor,

$$\begin{equation}{[{\mathrm{d}}^{2}S/\mathrm{d}{t}^{\prime 2}({t}_{s})]}^{-\eta /2-1/2}={[\mathbf{p}+\mathbf{A}({t}_{s})]}^{-\eta /2-1/2},\end{equation} \tag{ 16 }$$

where $\eta =Z/\sqrt{2{I}_{p}}$ is an effective principal quantum number of the initial state, cf reference [12]. The electron then follows a stationary trajectory specified by

$$\begin{equation}{\mathbf{r}}_{s}(\mathbf{p},t)=\mathbf{p}(t-{t}_{s})+{\Delta}\mathbf{G}(t,{t}_{s}),\end{equation} \tag{ 17 }$$

which has the following properties:

$$\begin{equation}{\mathbf{r}}_{s}(\mathbf{p},{t}_{s})=0,\quad {\dot {\mathbf{r}}}_{s}^{2}(\mathbf{p},{t}_{s})=-2{I}_{p},\quad {\dot {\mathbf{r}}}_{s}(\mathbf{p},{t}_{f})=\mathbf{p},\end{equation} \tag{ 18 }$$

where the first and second identities in equation (18) state that the electron starts at the origin with imaginary kinetic momentum at the complex interaction time, t_s , while the third identity states that the final kinetic and canonical momenta are equal at the end of the pulse.

Following Popruzhenko [12], we now apply a constraint on the momentum, p → p_nℓ,s , by imposing that the final squared angular momentum,

$$\begin{equation}{L}^{2}({t}_{f})={[{\mathbf{r}}_{s}(\mathbf{p},{t}_{f}){\times}\mathbf{p}]}^{2}={p}_{\perp }^{2}{\Delta}G{({t}_{f},{t}_{s})}^{2}=\ell (\ell +1),\end{equation} \tag{ 19 }$$

of the stationary trajectory has the quantized squared angular momentum, ℓ(ℓ + 1), of the target Rydberg state. Given that the angular momentum condition, equation (19), determines the square of the transverse momentum component,

$$\begin{equation}{p}_{\perp }^{2}=\frac{\ell (\ell +1)}{{\Delta}G{({t}_{f},{t}_{s})}^{2}},\end{equation} \tag{ 20 }$$

the component p_z parallel to the laser field can now be solved using equation (15) as

$$\begin{equation}{p}_{z}=-A({t}_{s}){\pm}i\sqrt{2{I}_{p}+{p}_{\perp }^{2}}=-A({t}_{s}){\pm}i\sqrt{2{I}_{p}+\frac{\ell (\ell +1)}{{\Delta}G{({t}_{f},{t}_{s})}^{2}}}.\end{equation} \tag{ 21 }$$

In this way the stationary momentum has been eliminated and one remaining condition is required to solve for t_s . Based on the LRL-vector, a condition on the final energy is imposed,

$$\begin{equation}\frac{{\mathbf{p}}^{2}}{2}-\frac{Z}{\sqrt{{\mathbf{r}}_{s}{(\mathbf{p},{t}_{f})}^{2}}}=-\frac{{Z}^{2}}{2{n}^{2}},\end{equation} \tag{ 22 }$$

that introduces a dependence on the principal quantum number, n, of the Rydberg state. The Rydberg matching conditions, equations (19) and (22), are used instead to the condition for stationary momentum, equation (13).

In order to find approximate solutions to the SFA-FT model, we reformulate equation (22) as

$$\begin{equation}A({t}_{s})={\pm}\sqrt{\frac{2Z}{\sqrt{{\mathbf{r}}_{s}^{2}}}-\frac{{Z}^{2}}{{n}^{2}}-\frac{\ell (\ell +1)}{{\Delta}{G}_{s}^{2}}}{\pm}i\sqrt{2{I}_{p}+\frac{\ell (\ell +1)}{{\Delta}{G}_{s}^{2}}},\end{equation} \tag{ 23 }$$

where p has been eliminated and r_s and ΔG_s are functions of t_s . Without loss of generality, we consider an electric field peak at t = 0, E(t) ≈ −E₀ cos ωt. Solutions to equation (23) are then expected to occur for t_s ≈ 0, such that A(t_s ) ≈ A₀ sin(ωt_s ) ≈ ωA₀ t_s . If the right-hand side of equation (23) is approximated with the zeroth-order interaction time, ${t}_{s}^{0}=0$, the first-order solutions become

$$\begin{equation}{t}_{s}=\frac{1}{\omega {A}_{0}}\left({\pm}\sqrt{\frac{2Z}{{r}_{0}}-\frac{{Z}^{2}}{{n}^{2}}-\frac{\ell (\ell +1)}{{\Delta}{G}_{0}^{2}}}{\pm}i\sqrt{2{I}_{p}+\frac{\ell (\ell +1)}{{\Delta}{G}_{0}^{2}}}\right).\end{equation} \tag{ 24 }$$

Here, the second term is always imaginary and it describes the tunneling process (the positive sign must be chosen for damping). In contrast, the first term describes the recapture process, and it may be either real or imaginary depending on the final radial position of the electron, ${r}_{0}=\sqrt{{\mathbf{r}}_{0}^{2}}$. Real contributions to t_s are found for

$$\begin{equation}{r}_{0}{< }\frac{2{n}^{2}}{Z+{n}^{2}\ell (\ell +1)/(Z{\Delta}{G}_{0}^{2})}{< }\frac{2{n}^{2}}{Z},\end{equation} \tag{ 25 }$$

which implies that the final electron position should not exceed a Kepler orbit with zero angular momentum. It is known from classical mechanics that increasing the angular momentum reduces the eccentricity of the orbit, and in the limit of large n and ℓ the orbit becomes circular with radius n²/Z. In the SFA-FT model we find that increasing the angular momentum, ℓ, indeed decreases the critical radius for recapture. This shows that the strict return condition for stationary momentum, given in equation (13), has been relaxed, by equations (19) and (22), to include returns to the entire volume of the final Rydberg state. Provided that this return condition is satisfied, there are two solutions for t_s , one on either side of the electric field maximum,

$$\begin{equation}\mathrm{R}\mathrm{e}({t}_{s})={\pm}\frac{\sqrt{2Z/{r}_{0}-{Z}^{2}/{n}^{2}-\ell (\ell +1)/{\Delta}{G}_{0}^{2}}}{\omega {A}_{0}},\end{equation} \tag{ 26 }$$

where the solution after (before) the electric field peak will be labeled α (β) in section 3. Further, it can be seen that an increased angular momentum, ℓ, decreases the real time delay Δt_s between the α and β solutions. At some critical angular momentum, ℓ > ℓ_n , the real contribution to t_s is prevented, leading to an imaginary contribution also by the first term on the right-hand side of equation (23). This signals that the electron trajectory ends up outside the recapture volume of the Rydberg state. In this case the positive sign must be selected also for the first term in equation (24), such that the stationary time leads to suppression of the recombination probability for that trajectory. The critical angular momentum can be estimated as

$$\begin{equation}{\ell }_{n}=-\frac{1}{2}+\sqrt{\frac{1}{4}+{\Delta}{G}_{0}^{2}\left(\frac{2Z}{{\Delta}{G}_{0}}-\frac{{Z}^{2}}{{n}^{2}}\right)},\end{equation} \tag{ 27 }$$

by setting the expression in the first square root of equation (23) to zero, and assuming that the motion is mostly due to the electric field so that r_s can be approximated by ΔG₀ = ΔG(t_f , 0). We note that this estimate gives a slightly smaller value for ℓ_∞ ≈ 6.84 compared to ℓ_∞ ≈ 7.31 as computed by equation (33) in reference [12] (for the pulse parameters considered in this article).

In general, there will be multiple electric field peaks in ultrashort laser pulses, which lead to multiple sets of complex interaction times, t_s that depend on the exact laser field. These solutions will be studied by solving equation (23) numerically, including the self-consistent dependence of t_s on the right-hand side, in the next section.

In this subsection we show solutions t_s to equation (23) that are obtained using Newton's method. Ultrashort laser pulses, with vector potential given in equation (3) for A₀ = 1.147 a.u. τ = 2T and period T = 2.67 fs, are considered. This corresponds to a full width at half maximum duration of 1.94 fs and 1.5 × 10¹⁴ W cm⁻² peak intensity.

3.1.1. Dependence on angular momentum of Rydberg state

In figure 1 we show the ℓ dependence of the stationary interaction time, t_s , for a hydrogen atom to a Rydberg state with principal quantum numbers n = 5 (a) n = 8 (b) and n = 20 (c), that is subject to an ultrashort laser pulse with cosine-like electric field (φ = 0). The curves for t_s are computed by continuous increase of ℓ from 0 to n − 1 with the quantized values, ℓ = 0, 1, ..., n − 1, marked by crosses. In figures 1(a)–(c) we observe four classes of solutions for t_s , labeled by α, β, γ and δ. The α and β solutions start far apart on the real axis with similar magnitude on the imaginary axis for ℓ = 0. The two solutions become closer together on the real axis as ℓ is increased. Clearly the α and β solutions correspond to the approximate solutions, t_s , for equation (23) found earlier in section 2.3. Since α starts after β, we may refer to the two classes of solutions as the short and long FT trajectories, respectively. It should be noted that these trajectories are not the same as the short and long HHG trajectories that both start after the electric field peak and both return to the origin for recombination with the atom. Because the short FT trajectory (α) is created after the electric field peak it is such a returning trajectory that may scatter with the atomic potential. The long FT trajectory (β) starts before the electric field peak, which implies that it has a net drift momentum away from the atom and never returns to the origin. Thus, the β solutions are most likely well described by the SFA model, while the α solutions may be strongly distorted due to rescattering with the ion. This point is illustrated in figure 2, where we show the trajectories associated to the different types of solutions displayed in figures 1(a) and (b), for n = 5, ℓ = 1 and n = 8, ℓ = 7. The trajectories are plotted from the exit point of tunneling, Re{t_s }, to the end of the pulse, t_f . In accordance with the discussion above, we see that in figure 2(a) the α trajectory crosses Re{r_z } = 0, while the β remains relatively far away from the core, but still within the recapture radius of ∼50 Bohr radii, estimated by equation (25) (the horizontal line in figures 2(c) and (d) indicate the tighter of the two bounds). In figures 2(b) and (d) on the other hand, we see that both α and β end up outside the recapture volume.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can be seen in figure 1(a) that the α and β solutions do not 'collide' in the complex plane for n = 5 as ℓ is increased from 0 to 4. This is because the angular momentum is restricted to low values below the critical angular momentum, ℓ₅ ≈ 4.51 > 4, estimated using equation (27) with ΔG₀ = 26.8 Bohr radii. In contrast, α–β collisions in the complex plane occur for both n = 8 (b) and n = 20 (c) somewhere in the range 6 < ℓ < 7, in good agreement with the estimates ℓ₈ ≈ 6.03 and ℓ₂₀ ≈ 6.84, respectively. Because the α solution has the correct signs for exponential damping when the electron is not recaptured classically, it is the physically correct solution for ℓ > ℓ_n . The β solution is an unphysical solution for ℓ > ℓ_n .

While the α and β solutions are consistent with the monochromatic field case, discussed in reference [12], the two additional solutions, γ and δ in figure 1, have not been reported previously. The γ and δ solutions have completely different behavior as compared with the α and β solutions. The real parts of γ and δ are of similar magnitude, while the imaginary parts of γ and δ differ in magnitude. The difference of imaginary parts tends to increase further with ℓ. Because the imaginary part of δ is smaller than that of α, β or γ it is essential to know if δ is a physically sound solution that should be included to model FT from ultrashort laser pulses. In order to answer this question, we consider the parallel canonical momentum, p_z (t_s ), shown in figures 3(a)–(c) with t_s taken from figures 1(a)–(c), respectively. We observe that both the α and β solutions have a mostly real p_z (t_s ) for low angular momentum, ℓ < ℓ_n . The α solution has a slightly larger magnitude of Re[p_z (t_s )], which implies a larger final speed of the shorter trajectory (α), when compared with the longer trajectory (β). Beyond, ℓ > ℓ_n , recapture is not possible and p_z , diverge in separate direction along the imaginary axis for the α and β solutions. In contrast, the γ and δ solutions have complex p_z (t_s ) for all possible angular momenta, which implies that they do not correspond to FT-like solutions, where the electron should evolve in the continuum before the recapture event. We do not have a physical understanding of the γ and δ solutions.

Zoom In Zoom Out Reset image size

In figures 1(d)–(f) we show the imaginary part of the semiclassical action of the solutions, Im[S(t_s )], computed with the t_s in ((a), (b) and (c)). This quantity leads to exponential damping of the processes. In the classically allowed region, ℓ < ℓ_n , we observe that the α and β solutions both have approximately the same weight. In the classically forbidden region, the weight of the α solution is reduced with increased angular momentum, ℓ > ℓ_n . In the case of n = 5 (g) and n = 8 (h), the weight of δ is larger than that of the γ weight and both α and β weights taken in the classical regime. This implies that inclusion of the δ solution would strongly alter the result of the SFA-FT model.

In order to ensure that only 'real' FT solutions are included in the model, we implement a multidimensional filter. This implies that FT trajectories that miss the recapture condition (α and β solutions with ℓ > ℓ_n ) will be removed, along with all non-FT solutions (γ, δ, etc solutions). The filter is defined as:

$$\begin{equation}\begin{aligned}\mathrm{I}\mathrm{m}\left\{{t}_{s}\right\}& { >}0\\ \mathrm{I}\mathrm{m}\left\{{S}_{0}({t}_{s})\right\}& {< }0\\ \vert \mathrm{I}\mathrm{m}\left\{{p}_{z}\right\}\vert & {< }0.1\vert \mathrm{R}\mathrm{e}\left\{{p}_{z}\right\}\vert \\ \vert {p}_{z}\vert & {< }0.5\quad \mathrm{a}.\mathrm{u}.,\end{aligned}\end{equation} \tag{ 28 }$$

where the first two conditions are standard for saddle point approximations in SFA, and they guarantee that | exp(−iS₀)| < 1. The third condition is specific to FT and it selects trajectories that have a classical recapture, while the final condition is used to remove unexpected solutions that originate for a small electric field, i.e. where E(Re{t_s }) is close to zero.

3.1.2. Dependence on carrier-envelope phase of laser pulse

The CEP of the laser pulse, φ, is defined in equation (12). Increasing the CEP leads to an advance of the electric field peak. In figure 4, it is shown how the CEP affects the solutions of equation (23), t_s (φ).

Zoom In Zoom Out Reset image size

Each individual panel shows t_s (φ) for a given ℓ, ranging from 0 to 7, and principal quantum number, n = 8. All obtained solutions are shown, including the classes: α, β, γ and δ, from the previous subsection. In order to identify the solutions that are real FT solutions, we mark solutions that pass the filter in false CEP color (see colormap). Solutions that are rejected by the filter are drawn with a consistent gray color.

First, we consider the simplest case, (n, ℓ) = (8, 0). At the centre of the pulse (t = 0 cycles), the four solutions: α, β, γ and δ have already been identified. The color scale indicates that for a given CEP, the short trajectory (α) is born later than the long trajectory (β). Physically, this corresponds to a laser pulse with an electric-field peak with a smaller phase delay. The α and β solutions can be followed back to the beginning of the pulse (t = −1 cycles) by increasing the CEP by more than π, thus, smoothly connecting classes of solutions in different half cycles of the laser pulse. Decreasing the CEP, which is the same thing as delaying the field peaks, shows that the short trajectory (α) stops being a real FT solution rather early in the pulse for t > 0.4 cycles. We observe that the long trajectory (β) remains a real FT solution for an extended time window in the pulse. It develops into non-real FT solutions for t > 0.8 cycles. In contrast, both γ and δ solutions are non-real FT solutions for all CEP in the entire time window of the pulse.

When angular momentum is increased, the picture becomes more complicated because additional solutions to equation (23) are found. In the case of (n, ℓ) = (8, 1) the decreased recapture volume, see equation (25), leads to α and β solutions becoming non-real FT solutions at the start of the pulse. At the centre of pulse both α and β are real FT solutions. The start of real FT trajectories on the α solutions is later than that of the β solutions. This indicates that the α trajectories are more sensitive to laser field envelope effects. In addition to the non-real γ and δ solutions, additional solutions are found at the start and end of the laser pulse, as well as at small imaginary time. All of these additional solutions are rejected by the filter as non-real FT solutions and consistently in gray color. Further, it can be seen that the filter selects the 'real' FT solutions, that is to say before the α and β solutions split up on the imaginary axis. With (n, ℓ) = (8, 5) we can clearly see the transition from real FT trajectories to the complex case, with the short trajectory (α) being the physical solution that decrease exponentially (by increase of its imaginary part), while the long trajectory can be traced to unphysically increased contribution. As the angular momentum is increased further, we observe that the time window for real FT trajectories of the α and β solutions, is reduced. While a small time window for real FT trajectories remains with (n, ℓ) = 8, 6, the highest angular momentum case, (n, ℓ) = (8, 7), shows that no real FT trajectories are possible for any CEP with the given ultrashort pulse.

To test if the SFA-FT theory is able to capture the physics that leads to strong field excitation, we present a comparison with results obtained by numerically solving the TDSE. The numerical propagation is done in a spectral basis, with the radial wavefunctions and matrix elements computed by B-splines, and the time propagation is done using a second-order differencing scheme [18] in the spectral domain. Propagation has been performed in both length and velocity gauge to ensure that the results are converged and gauge invariant. Figure 5 shows a CEP scan of n-populations that was made using the SFA-FT theory and the TDSE. To remove unphysical contributions in the SFA calculations, we use the filter from equation (28). The left and middle column were computed with the SFA theory, including quantum diffusion, using only the angular part ${Y}_{\ell m}^{{\ast}}(\mathbf{p}/\sqrt{{p}^{2}})$ (left) or the full ${~{\phi }}_{n\ell m}^{{\ast}}(\mathbf{p})=\left\langle {\phi }_{n\ell m}\vert \mathbf{p}\right\rangle$ (middle) for the projection prefactor in equation (10). The right column shows the results of the TDSE simulations. The top, middle and bottom rows were calculated for two, three and four cycle pulses, i.e. τ = 2T, 3T and 4T in equation (3), respectively. The populations have been normalized to the maximum of each subplot, so that the same color scale can be used for all plots. The sharp steps that are seen in the SFA results of figure 5, are due to the filtration process. As shown in figure 4, the evolution of the SFA-FT solutions is very complicated, which makes filtration difficult.

Zoom In Zoom Out Reset image size

As the CEP is varied we see an anti-modulating structure in the TDSE results where a population maximum for one n = 5 value is accompanied by minima for other n = 6, and vice versa. This can be clearly seen in figure 5(c) at φ = 0.5. This behavior is not found in the SFA-FT results. When n ⩾ 6, the modulations for different n have similar shapes in both TDSE and SFA. In general the qualitative agreement between SFA-FT and TDSE is better if we use the momentum-space Rydberg state as the projection prefactor, ${~{\phi }}_{n\ell m}^{{\ast}}(\mathbf{p})$. For instance, the full wavefunction is needed to find the correct decrease in population as n increases.

Figure 6 shows the results for n = 5 from figure 5, resolved in ℓ and summed over all ℓ. A qualitative agreement between the total population for n = 5 is observed between the model ((b), (e) and (h)) and the TDSE ((c), (f), (i)). A maximum is found for a cosine-like pulse (close to φ = 0) with two and four cycles, while a minimum is found for three cycles. The strength of the peaks as a function of CEP are not correctly described by the model, and it is found that the projection factor is essential for the qualitative agreement with TDSE, see figures 6(a), (d) and (g). The fact that we see several maxima as a function of CEP, is a sign of quantum interference between different trajectories. In contrast, if the contributions from different trajectories are summed incoherently a single smooth maximum is observed in the given CEP range, as expected from the tunneling step of the process (not shown). From these plots it is clear that although we can find qualitative agreement between the SFA-FT and TDSE results when the results are summed over ℓ, the contributions from individual ℓ channels are completely different in the two approaches. In the TDSE results, the total population for n = 5 comes mostly from ℓ ⩾ 2, with a tendency towards being completely dominated by ℓ = 4 for the longer pulses. For both SFA calculations we mostly see contributions from ℓ ⩽ 2, with ℓ = 4 barely contributing. There is also a difference between the two SFA prefactors. When we only use the angular part of the momentum space wavefunction, the dominant contributions come from ℓ = 1 and ℓ = 2, while when we include the full wavefunction the dominant contributions come from ℓ = 0 and ℓ = 1. This can explain why we see larger discontinuities for the total population in n = 5 when we only include the angular factor, since more solutions are filtered out as ℓ increases, see figure 4. Excitation to n = 5 corresponds to 8.42 laser photons, which is in between even and odd photon numbers. With ultrashort pulses the width of the photons allows for overlapping parity to the final state. For the shortest pulses considered, even and odd parity states have comparable populations, see figure 6(c), while for the longest pulses considered, even parity is favored, see figure 6(i). Similar mixing of parity is observed also in the SFA-FT model.

Zoom In Zoom Out Reset image size

In order to better understand the role that the continuum states play in the strong-field population of excited states, we performed TDSE simulations where we either remove them entirely from the state space, or add an imaginary part to the energy of each continuum state so that their amplitudes become damped over a subcycle timescale during propagation. These calculations were performed in length gauge. We found that the continuum states play a crucial role in getting the correct excited state distribution, and it appears to be especially important for n ⩾ 4. Figure 7 shows a comparison between the results of standard propagation, removing the continuum states, and changing the energy of the continuum states to E → E − 0.5i. The results for n = 2 and n = 3 have not been affected very much by the addition of complex continuum energies, whereas the higher excited states show a significant difference to the results of the standard propagation. Propagating without the continuum states predicts a maximum at n = 3 instead of n = 5, but the yield to higher n is similar. This suggests that the FT model can have qualitative agreement with the TDSE results for n ⩾ 4, and why the same is not true for n = 2, 3. We find similar results when the the imaginary part of the energy is 1 and 0.1.

Zoom In Zoom Out Reset image size