A Three-step Model of High Harmonic Generation using Complex Classical Trajectories

TL;DR

This paper introduces a novel complex trajectory approach to high harmonic generation that unifies classically allowed and forbidden processes, providing a comprehensive dynamical framework for understanding strong field physics effects.

Contribution

It presents a new trajectory formulation that models both allowed and forbidden processes in high harmonic generation within a single, seamless dynamical system.

Findings

Complex trajectories orbit the nucleus, forming the Coulomb ground state.

The model describes tunnel ionization, unbound evolution, and recollision seamlessly.

The approach offers mechanistic insights into strong field physics phenomena.

Abstract

We present a new trajectory formulation of high harmonic generation that treats classically allowed and classically forbidden processes within a single dynamical framework. Complex trajectories orbit the nucleus, producing the stationary Coulomb ground state. When the field is turned on, these complex trajectories continue their motion in the field-dressed Coulomb potential and therefore tunnel ionization, unbound evolution and recollision are described within a single, seamless framework. The new formulation can bring mechanistic understanding to a broad range of strong field physics effects.

Equations

$${\Psi }_0(x)={\mathrm{e}}^{\frac{\mathrm{i}}{\hslash }{S}_0(x)}{S}_0(x),x\in \mathbb{C}.$$

$$p_0(q_0)={\frac{\partial S_0(x)}{\partial x}|}_{x=q_0},$$

$$\begin{array}{rl}\dot{q}& =\frac{1}{m}p\\ \dot{p}& =-V^{\prime }(q)\\ \dot{S}& =\frac{1}{2m}{p}^2-V(q)+\frac{\mathrm{i}\hslash }{2m}\alpha \\ \dot{\alpha }& =-V^{\prime \prime }(q)-\frac{1}{m}{\alpha }^2\end{array}$$

$$2\gamma q_t(q_0)-\mathrm{i}{p}_t(q_0)$$

$$\sigma (q_0)$$

$$⟨x|g_{\gamma }(\xi (q_0))⟩={\left(\frac{2\gamma }{\pi }\right)}^{\frac{1}{4}}{\mathrm{e}}^{-\gamma {\left[x-{\bar{q}}_t(q_0)\right]}^2-\mathrm{i}{\bar{p}}_t(q_0)\left[x-{\bar{q}}_t(q_0)\right]}.$$

$${\Psi }_t(x)=\int \frac{-\mathrm{d}{q}_0}{4\pi \gamma }|J(q_0)|⟨x|g_{\gamma }(\xi (q_0))⟩⟨g_{\gamma }(\xi (q_0))|\mathcal{U}|{\Psi }_0⟩.$$

$$\frac{\mathrm{d}\xi (q_0)}{\mathrm{d}{q}_0}=\left(2\gamma M_{qp}-\mathrm{i}{M}_{pp}\right){\alpha }_0(q_0)+\left(2\gamma M_{qq}-\mathrm{i}{M}_{pq}\right).$$

$$a(t)=⟨\Psi (t)|\ddot{x}|\Psi (t)⟩=\int \mathrm{d}x|\Psi (x,t){|}^2\frac{\partial V(x,t)}{\partial x}.$$

$${\dot{M}}_{\mathrm{pp}}$$

$${\dot{M}}_{\mathrm{qp}}$$

Full text

strong field Coulomb dynamics, high harmonics spectra, semiclassics, complex trajectories

Introduction

The strong field dynamics of Coulombic systems has been at the center of attention of the atomic, molecular and optical physics communities for decades. Particularly fascinating are the effects of above threshold ionization (ATI) Agostini et al. (1979) and high harmonic generation (HHG) McPherson et al. (1987); Ferray et al. (1988); Li et al. (1989). Besides being of interest for its fundamental physics, HHG is a powerful tool for generating XUV pulses that can be used to further investigate atomic and molecular quantum dynamics.

Above threshold ionization and high harmonic generation are intrinsically quantum mechanical phenomena, and where numerical evaluation of the time dependent Schrödinger Equation (TDSE) is possible it faithfully reproduces the experimentally observed quantities. However, numerical simulation does not always provide mechanistic understanding, and hence there is great interest in trajectory-based classical and semi-classical models to describe these phenomena Protopapas et al. (1996); van de Sand and Rost (1999).

Consider the three-step model of HHG introduced by Corkum Corkum (1993) (see also Kulander et al. (1993)). In the three-step model, the HHG process is decomposed as follows: 1) tunnel ionization 2) free propagation in the presence of the field, and 3) recollision and recombination producing HHG emission Corkum (1993); Lewenstein et al. (1994). It is found that classical models are accurate for the second step of the three-step model – the free propagation in the presence of the field. But real-valued trajectories cannot describe the tunnel ionization that takes place in the first step. Moreover, real-valued trajectories in the Coulomb system cannot describe a stationary initial wavefunction because all trajectories fall into the nucleus. This undermines a consistent formulation for sampling initial conditions from the stationary wavefunction 111This is true for 1D models and for multi electron systems. Two and three dimensional single electron models do not have this problem.. As a result of these limitations, classical-based models have generally chosen initial conditions on the outside of the field-induced Coulomb barrier, with the choice of initial momentum zero.

In this paper we show that complex-valued classical trajectories can overcome all of the above problems, accounting for tunnel ionization, laser dressed electron propagation and recollision, all within a single consistent framework. For the purposes of presentation we define a modified three-step model:

1. Construction of the stationary ground state Coulomb wavefunction using complex trajectories. It is difficult to overstate how non-trivial this step is, as will become clear below;
2. Evolution of the complex trajectories in the presence of the Coulomb potential plus time-dependent strong field. Note that in this step tunneling ionization, free evolution and recombination are all treated in a seamless, unified way.
3. Construction of the high harmonic spectrum from the Fourier transform of the time-dependent dipole acceleration. At every stage, the time-dependent wavefunction is reconstructed and shown to be in semi-quantitative agreement with the exact quantum result.

This is by no means the first attempt to use complex trajectories to describe strong field ionization Perelomov et al. (1966, 1967); Perelomov and Popov (1967); Keldysh (1964); Faisal (1973); Reiss (1980); Salières et al. (2001); Shafir et al. (2012); Kaushal and Smirnova (2013); Torlina et al. (2015); Pisanty and Ivanov (2016); Kaushal and Smirnova (2018) but we believe this is the first treatment where all portions of the three-step model are described within a single, consistent trajectory framework, sampled from an initial Coulomb eigenstate. In some sense it is the completion of the program of the classic papers by Perelomov, Popov, and Terent’ev (PPT) Perelomov et al. (1966), Corkum Corkum (1993) and Lewenstein et al. Lewenstein et al. (1994).

Theory

The initial state of the system $\Psi_{0}(x)$, is analytically continued into the complex plane. This defines a complex-valued action $S_{0}(x)$

[TABLE]

A manifold of complex-valued initial conditions $q_{0},p_{0}\in{\mathbb{C}}$ is defined through the relation

[TABLE]

with a stability parameter $\alpha_{0}(q_{0})=\left.\frac{\partial^{2}S_{0}(x)}{\partial x^{2}}\right|_{x=q_{0}}$. The complex trajectories are then propagated by integrating the equations of motion (EOM)

[TABLE]

where $m$ is the mass of the system. The complex action $S=S^{\rm cl}+S^{\rm qm}$, where $S^{\rm cl}=\int\mathrm{d}t\frac{p^{2}}{2m}-V(q)$ and $S^{\rm qm}=\int\mathrm{d}t\frac{{\mathrm{i}}\hbar}{2m}\alpha$.

In the final value coherent state propagator (FINCO) method Zamstein and Tannor (2014), each trajectory is associated with a matrix element of the propagator in the coherent state representation $\left<g_{\gamma}(\xi(q_{0}))\vphantom{\mathcal{U}\Psi_{0}}\right|\mathcal{U}\left|\Psi_{0}\vphantom{g_{\gamma}(\xi(q_{0}))\mathcal{U}}\right>$. Equations  (3) subject to initial conditions (2) determine the complex coordinate $q_{t}(q_{0})$ and momentum $p_{t}(q_{0})$ at time $t$. This in turn provides a coherent state label $\xi(q_{0})$ for the bra Gaussian $\left\langle g_{\gamma}(\xi(q_{0}))\right|$ through the relation

[TABLE]

where $q_{t},p_{t}\in{\mathbb{C}}$, $\bar{q}_{t},\bar{p}_{t},\gamma\in{\mathbb{R}}$ Huber and Heller (1987). The coherent state matrix element is given by $\left<g_{\gamma}(\xi(q_{0}))\vphantom{\mathcal{U}\Psi_{0}}\right|\mathcal{U}\left|\Psi_{0}\vphantom{g_{\gamma}(\xi(q_{0}))\mathcal{U}}\right>=\phi(q_{0})\,\mathrm{e}^{\sigma(q_{0})}$, where the exponent $\sigma(q_{0})$ has the form

[TABLE]

Associated with this bra Gaussian is a ket Gaussian centered at $\xi^{\ast}=2\gamma\bar{q}_{t}(q_{0})+{\mathrm{i}}\bar{p}_{t}(q_{0})$,

[TABLE]

The time evolved wavefunction $\Psi_{t}(x)$ is represented in a basis of these ket Gaussians

[TABLE]

The Jacobian $|J(q_{0})|=\left|\frac{\mathrm{d}\xi(q_{0})}{\mathrm{d}q_{0}}\right|^{2}$ and prefactor $\phi(q_{0})=\left(8\gamma\pi\right)^{\frac{1}{4}}\left[\frac{\mathrm{d}\xi(q_{0})}{\mathrm{d}q_{0}}\right]^{-\frac{1}{2}}$, where

[TABLE]

The $M_{ab}$ are stability matrix elements defined by $M_{ab}=\frac{\partial a_{t}(q_{0})}{\partial b_{0}(q_{0})}$ with $a,b\in\{p,q\}$. EOM for the $M_{ab}$ are given in App. A.

Complex time

When the EOM of the previous section are applied to the Coulomb system, the methodology fails at the first step: no meaningful reconstruction of the ground state eigenfunction is obtained. Real-valued bound state initial conditions all fall into the Coulomb singularity. Strictly complex-valued trajectories, on the other hand, diverge, as the complex Coulomb potential is repulsive everywhere except on the real axis. Hence a stable ground state can not be constructed from either real or complex trajectories.

The resolution to this conundrum comes from allowing time to be complex. The complex time plane has singularities – points where the classical momentum goes to infinity in a finite time. Paradoxically, these singularities, rather than being a nuisance, provide the solution to the problem by changing the topology of the trajectories. For instance, it is known that circumnavigating a time-singularity in barrier scattering problems can transform a transmitted trajectory into a reflected trajectory Doll et al. (1973); Kay (2013); Petersen and Kay (2014). In the Coulomb system, circumnavigating complex time singularities leads to complex trajectories that orbit the nucleus. Technically, these singularities cause the solution $q_{t}(q_{0}),p_{t}(q_{0})$ to be multivalued, but the motion is continuous when viewed as evolving on a multi-sheeted Riemann surface. Motion on these sheets may correspond to either classically allowed or classically forbidden processes, depending on the trajectory initial conditions.

To understand more fully the influence of circumnavigating a singularity in time, we note that the singularities in complex time arise from singularities in complex position Doll et al. (1973); Kay (2013). In Ref. Koch and Tannor (2018a) a full calculus was introduced relating the asymptotic approach towards a singularity of the potential to the singularities in complex time. As shown in Koch and Tannor (2018a), the Coulomb singularity at the origin gives rise to three Riemann sheets for each singularity encounter in time. Controlled circumnavigation of the singularity in time leads to complex trajectories that cycle through these three Riemann sheets, corresponding to trajectories that orbit the nucleus (Fig. 1).

Associated with the singularities in complex position there are always caustics, although the mathematical connection between these two is yet to be discovered. These caustics lead to divergence of a portion of the complex trajectories, those that lie beyond what are known as Stokes lines. A rigorous procedure for removing these unphysical trajectories was introduced in Ref. Koch and Tannor (2018b).

The three-step model. Step 1: Constructing the ground state of the Coulomb system

The first step is to describe the stationary initial state of the Coulomb system in terms of orbiting trajectories. Consider the anti-symmetric, Cartesian 1D Coulomb ground state $\Psi_{0}(x)=2x\,\mathrm{e}^{-|x|}$. The strategy is to let these trajectories orbit the nucleus freely until the field is turned on and then seamlessly continue their propagation in the presence of the field.

A complex time contour that loops around a singularity time passes onto a second Riemann sheet where the corresponding trajectory orbits the nucleus. On this alternate sheet, an additional singularity time is present which is in turn circumnavigated by a loop contour. Continuing in this way, each loop corresponds to a trajectory passing from one sheet of the Riemann surface to another, followed by propagation on that new sheet until the next singularity encounter. Figure 1 shows this transition between sheets for two successive singularity times.

The total complex time contour thus consists of a series of loops around successive singularity times as shown in Fig. 2.

Note that the wavefunction is always reconstructed at purely real times. The necessary excursions back to the real time axis are not shown in Fig. 2.

The locations of singularity times depend on the trajectory initial conditions. Propagating the entire manifold with suitably chosen contours yields the desired stationary ground state depicted in Fig. 3, the only time-dependence being a coordinate-independent phase rotation.

Both the magnitude and phase of the wave function are in semiquantitative agreement with the exact quantum results.

Step 2: Strong field dynamics of the Coulomb system

Having described the construction of the initial wavefunction, we now introduce the time-dependent electric field. The laser field is modeled in the dipole approximation and length gauge by the additional, time dependent potential $\tilde{V}(q,t)=eqF_{0}\sin(\omega t)$ with the field strength $F_{0}=7.35\cdot 10^{-2}\,\mathrm{au}$, electronic charge $e$ and frequency $\omega=7.35\cdot 10^{-2}\,\mathrm{au}$. The Keldysh parameter is $\gamma=1$.

The integration contour in the presence of the field is shown in Fig. 4. Note that the addition of $\tilde{V}(q,t)$ does not introduce additional singularities in complex space. It does, however, shift the singularity times and introduces additional encounters with the Coulomb singularity. The loops around the shifted singularity times correspond to orbiting trajectories that remain near the core of the nucleus, modified by the presence of the field 222In practice, in order to avoid numerical inaccuracies in the ground state, we skip the initial field free loops and start with the field dressed loops at $t=0$.. However, the loops that circumnavigate the new singularity encounter correspond to a new kind of trajectory dynamics: ionization of the electron and subsequent recollision with the nucleus.

The ionizing and recolliding trajectories result from three distinct parts of a complex time contour :

1. A series of singularity loops causing the trajectory to orbit the singularity at the origin. Termination of the loop sequence initiates tunnel ionization.

2. Straight line integration in complex time results in a field-mediated excursion from the nucleus. Depending on the trajectory initial conditions, the trajectory will ionize or return to the core at a later time.

3. A loop around a new singularity time that is not present in the field free case induces recollision with the nucleus.

Each additional cycle of the field results in additional contributions to the emitted wavepacket. These contributions are generated via contours analogous to those of the first ionization cycle, albeit with a varying number of initial loops. E.g. the first ionization cycle begins after three initial loops whereas the second ionization cycle requires fifteen loops. Note that the time of maximum field strength roughly coincided with three periods of phase rotation of the ground state and a full cycle of the field corresponds roughly to twelve periods of rotation of the ground state such that the second field maximum occurs after about fifteen periods of the ground state rotation. Terminating the initial loop series earlier or later yields trajectories that describe the same process but that disagree with the quantum results. The physical significance of these alternate contours is not yet understood.

Adding the contributions from all distinct contours that make significant numerical contributions yields the wavefunction in Fig. 5.

The ionized part of the wave function plotted for 1.5 periods of the field shows excellent agreement with the quantum result, including the complex nodal pattern resulting from the superposition of the three distinct contributions, one from each of the three half cycles of the field prior to the evaluation time.

As the laser field changes sign, the previously emitted part of the wavefunction is partially redirected towards the nucleus. This recolliding wavefunction recombines with the ground state of the system, and can be computed from the trajectory based reconstruction.

Step 3: Calculating the high harmonic spectrum

The spectrum of the radiation emitted as the recolliding wavefunction recombines with the ground state of the system can be calculated as the Fourier transform of the dipole acceleration computed in the form Burnett et al. (1992)

[TABLE]

A comparison of the quantum and complex trajectory calculation is shown in Fig. 6. The locations of the peaks as well as the plateau structure, including the cutoff, are reproduced very well. Remaining deviations in the intensity of the peaks are due to errors in the reproduction of the core wavepacket.

Several comments about our numerical procedure. The FINCO wavefunction was reconstructed from four distinct classical processes: the first and third half cycle ionization and the second and fourth half cycle recollision. The dominant regions in the initial manifold contributing significantly to each of these processes were selected manually. Since the quantitative simulation of the core wavepacket is somewhat more involved, we added to the ionized and recolliding parts of the wavefunction, the Coulomb ground state with its time-dependent phase, $\Psi(x,t)\rightarrow\Psi(x,t)+\Psi_{0}(x)\,\mathrm{e}^{\frac{{\mathrm{i}}}{2}t}$.

All of the data presented was computed from initial conditions in only the positive half of coordinate space, evolving on a single Riemann sheet of the potential $V(q)=\mp\frac{1}{q}$. Instead of computing the negative half from the initially antisymmetric ground state $\Psi_{0}(x)$ evolving on the other sheet of the potential, the dipole acceleration was post-symmetrized according to $a(t)\rightarrow a(t)-a(t-\frac{\pi}{\omega})$.

Summary

To summarize, we have shown how strong field dynamics of the Coulomb system, including high harmonic generation, can be modeled using complex classical trajectories. In this approach the spectrum emerges in the three following steps:

1. Trajectories orbit the singularity in the complex plane, generating the Coulomb ground state.
2. Under the influence of the strong field, a subset of the trajectory manifold ionizes, with a subset of the subset recolliding with the nucleus.
3. Recombination with the nuclear wavefunction yields high harmonic radiation. There is no need for a separate tunnel ionization step as in the traditional three step model; it is included automatically in the second step.

The high harmonic spectrum reconstructed from the FINCO trajectories agrees well with the full quantum results. With the complex trajectories correctly describing tunnel ionization of the Coulomb ground state, a number of exciting future possibilities present themselves.

There has been significant debate in the literature whether the transition of the tunneling electrons through the barrier happens instantaneously or with a delay Kaushal and Smirnova (2013); Landsman and Keller (2014); Zipp et al. (2014); Pedatzur et al. (2015); Torlina et al. (2015); Ni et al. (2016). Since the complex trajectories used in the present work are rigorously derived from the Coulomb ground state without artificially imposing tunnel exit initial conditions, tunneling times may be read off directly. Note that each trajectory has its own independent tunneling time, allowing for transparent analysis and interpretation.

Aside from the manual selection of the significant regions of the initial manifold, the only approximation in the present approach is a truncation of an infinite hierarchy of EOMs at second order Goldfarb et al. (2006), corresponding to a saddle point evaluation of the path integral Schiff et al. (2011). Since no assumption is made on physical parameters such as the strength of the field, we expect that the complex trajectory method presented here will provide mechanistic insight into a broad range of strong field effects in atomic and molecular systems.

Acknowledgements We thank the organizers and participants of the Advanced Study Group “Semiclassical Methods: Insight and Practice in Many Dimensions” and the Focus Workshop “Quantum and Semiclassical Trajectories”, both held at the MPIPKS of Dresden, for many stimulating discussions. Financial support from the Israel Science Foundation (1094/16) and the German-Israeli Foundation for Scientific Research and Development (GIF) are gratefully acknowledged.

Appendix A Equations of motion for stability matrix elements

The stability matrix elements $M_{ab}=\frac{\partial a_{t}(q_{0})}{\partial b_{0}(q_{0})}$ with $a,b\in\{p,q\}$ obey the following EOMs

[TABLE]