Communication: Phase space approach to laser-driven electronic wavepacket propagation

We propose a phase space method to propagate a quantum wavepacket driven by a strong external field. The method employs the periodic von Neumann basis with biorthogonal exchange recently introduced for the calculation of the energy eigenstates of time-independent quantum systems [A. Shimshovitz and D. J. Tannor, Phys. Rev. Lett. (in press) [e-print arXiv:1201.2299v1]]. While the individual elements in this basis set are time-independent, a small subset is chosen in a time-dependent manner to adapt to the evolution of the wavepacket in phase space. We demonstrate the accuracy and efficiency of the present propagation method by calculating the electronic wavepacket in a one-dimensional soft-core atom interacting with a superposition of an intense, few-cycle, near-infrared laser pulse and an attosecond extreme-ultraviolet laser pulse.

With the emergence of attosecond laser technology,¹ there is the fascinating prospect of observing and controlling the correlated dynamics of multiple electrons on its natural time scale of ten to one hundred attoseconds.² In order to unravel the complex and sometimes counter-intuitive^3,4 quantum dynamics from the experimental data, and to develop theories that reproduce the essence of the dynamics,⁵ an accurate and efficient numerical method to simulate the multi-electron wavepacket dynamics is indispensable.

However, accurate simulation of the electronic dynamics in a high-intensity laser field is a challenging task: the electronic wavepacket is dispersed by the laser field over a wide region of coordinate space while retaining high momentum near the atomic nuclei. Straightforward representation of the wavefunction on a equally-spaced coordinate grid [i.e., a Fourier grid (FG)] requires a large range with a small interval between points. Simulation on such a large grid quickly becomes prohibitive as the number of degrees of freedom (DOF) increases. Even with sophisticated treatments such as multi-configuration time-dependent Hartree-Fock (MCTDHF),⁶ simulation of ionization dynamics have been limited to small systems such as the helium atom and the hydrogen molecule.^7–12 

In this article, we present a new approach to solving the time-dependent Schrödinger equation (TDSE), based on a phase space perspective. The resultant propagation method is simple, accurate, stable, and efficient. As a first demonstration, we simulate the electronic wavepacket of a one-dimensional (1D) model atom in the combined fields of a high-intensity near-infrared (NIR) laser pulse and an attosecond extreme-ultraviolet (XUV) laser pulse.

In our approach, we utilize the localized nature of phase space Gaussians to prune the basis. However, a basis set of phase-space localized states is perforce non-orthogonal,¹³ and this has created difficulties in previous attempts^14–19 to truncate the basis to cover only the phase space region of dynamics. This long standing problem was solved recently, by the so-called periodic von Neumann (pvN) basis.^20,21 The pvN basis is generated from a set of Gaussians whose centers form a finite lattice in phase space. By imposing periodic boundary conditions on these Gaussians, the pvN basis becomes formally equivalent to the simple and accurate FG representation.²² Then, the basis biorthogonal to the pvN basis (henceforth the ‘pvb’ basis) can be used for a compact representation of a quantum state. This framework was successfully applied to the calculation of quantum energy eigenstates.²¹ Here, the framework is extended for solving the TDSE. Our strategy is to keep the individual elements in the basis set time-independent but to truncate the pvb basis in a time-dependent manner. This avoids the problem of a moving basis¹⁵ that can become over-complete as time elapses and can become unstable. In contrast, here we obtain a stable set of linear ordinary differential equations (ODEs) for the expansion coefficients in the truncated pvb basis.

Before explaining the pvN and pvb bases, we review the formalism of the FG basis to establish notation. We write the Fourier pseudo-spectral and spectral bases^23–25 as {|θ_m⟩}_{m = 1, …, N} and {|ϕ_m⟩}_{m = 1, …, N}, respectively. These bases are orthonormal and span the same N-D Hilbert space denoted here as

$H$⁠, resolving the identity in

$\mathcal {H}$

$H$ as

The bases {|θ_m⟩} and {|ϕ_m⟩} are localized at the FG points {x_m}_{m = 1, …, N} in the position space and {p_m}_{m = 1, …, N} in the momentum space, respectively. For any quantum state

$|Ψ⟩∈H$⁠,

$\langle \theta _m|\Psi \rangle = \langle x_m|\Psi \rangle \sqrt{\Delta x}$

$⟨θm|Ψ⟩=⟨xm|Ψ⟩Δx$ and

$\langle \phi _m|\Psi \rangle\break = \langle p_m|\Psi \rangle \sqrt{\Delta p}$

$⟨ϕm|Ψ⟩=⟨pm|Ψ⟩Δp$⁠, where Δx and Δp are the grid intervals in position and momentum spaces, respectively.

The pvN basis^20,21 

${|g̃j⟩}j=1,⋯,N$ is defined as

where {|g_j⟩}_{j = 1, …, N} are the phase space Gaussians,

whose centers {(q_j, p_j)}_{j = 1, …, N} constitute a finite lattice in the phase space with the unit cell of area 2πℏ and the momentum-to-position aspect ratio ℏγ. Note that the number of Gaussians N is the same as the number of FG points used.

The pvb basis²¹ 

${|b̃j⟩}j=1,⋯,N$ is defined to be biorthogonal (dual) to the pvN basis, i.e.,

$\langle \tilde{b}_l|\tilde{g}_j\rangle = \delta _{lj}$

$⟨b̃l|g̃j⟩=δlj$⁠. This gives

where S⁻¹ is the inverse of the overlap matrix

$Slj=⟨g̃l|g̃j⟩=Δx∑m=1N⟨gl|xm⟩⟨xm|gj⟩$ of the pvN basis. The pvN and pvb bases span the same Hilbert space

$\mathcal {H}$

$H$ as {|ϕ_m⟩} and {|θ_m⟩}, resolving the identity as

Thus, any wavepacket

$|Ψ(t)⟩∈H$ can be represented as

Due to the localized nature of the Gaussians {|g_j⟩}, the magnitude of

$⟨g̃j|Ψ(t)⟩$ can be extremely small if the corresponding classical system can not reach the phase space region around (q_j, p_j). Defining a set

$\mathcal {A}$

$A$ such that

$|\langle \tilde{g}_j|\Psi (t)\rangle |$

$|⟨g̃j|Ψ(t)⟩|$ is negligible if

$j\not\in \mathcal {A}$

$j∉A$⁠, we can approximate the wavepacket by a subset of the pvb basis as

where

$cj(t):=⟨g̃j|Ψ(t)⟩$⁠. Note that the set

$\mathcal {A}$

$A$ of active indices can be changed in time in order to keep the number

$N_{\mathcal {A}}$

$NA$ of elements in

$\mathcal {A}$

$A$ small at all time.

By substituting Eq. (7) to the TDSE, we obtain a set of linear ODEs for the active pvb coefficients

${cj}j∈A$⁠,

where Ω⁻¹ is the inverse of the overlap matrix

$Ωjl=⟨b̃j|b̃l⟩=(S−1)jl$ of the active pvb basis, and H(t) is the Hamiltonian operator of the system. The overlap and Hamiltonian matrix elements in Eq. (8) can be computed simply via the representations in {|θ_m⟩} and {|ϕ_m⟩}.^21,23 The matrix Ω⁻¹ is Hermitian positive-definite, and the elements

$\lbrace \langle \tilde{b}_l| H(t) |\tilde{b}_m\rangle \rbrace$

${⟨b̃l|H(t)|b̃m⟩}$ constitute an Hermitian matrix. Therefore, the product of these matrices yields all real eigenvalues,²⁶ and Eq. (8) can be solved stably by many standard numerical algorithms.

To demonstrate the accuracy and efficiency of the present method, we solve Eq. (8) for the electronic wavepacket of a 1D atom in the combined field of NIR and XUV laser pulses. The Hamiltonian of this system is given as

where H₀ is the field-free Hamiltonian expressed as

Here μ = 1 a.u. is the electron mass, e = −1 a.u. is the electron charge, −Qe = 1 a.u. is the charge of the atomic nucleus, a = 1 a.u. is the soft-core parameter, and ε₀ = 1/4π a.u. is the electric constant. The laser-electron coupling V(t) is, in the velocity gauge,

where A_NIR(t) and A_XUV(t) are the vector potentials of the NIR and XUV laser pulses, respectively. We used an NIR pulse of wavelength 800 nm, peak intensity 5 × 10¹³ W/cm², and duration 1.5 cycles (4 fs, FWHM of intensity profile). The XUV pulse had wavelength 15 nm, peak intensity 1 × 10¹² W/cm², and duration 5.0 cycles (250 as). The peak of the XUV pulse was delayed from that of the NIR pulse by 0.25 NIR cycles, as shown in Fig. 1.

The initial state was chosen as the ground state of H₀ with energy eigenvalue −0.66978 a.u., and the wavepacket was propagated from the turn-on of the NIR laser pulse at t_min to its end at t_max by the short-iterative Arnoldi algorithm^27,28 in a 6D Krylov space with a constant time-step of Δt = 0.0379 a.u. We divided the time span from t_min to t_max into 8 time segments and changed the active set

$A$ from one segment to the next. The number of FG points was N = 4096, and they were distributed over −750 a.u. ⩽ x ⩽ 750 a.u. and −8.58 a.u. ⩽ p ⩽ 8.58 a.u. This phase space rectangle was divided into 64 × 64 cells of aspect ratio ℏγ = 1.14 × 10⁻² a.u., each of which contained a Gaussian center (q_j, p_j).

In Fig. 2, snapshots of |c_j(t)|² are shown by the color scale of the ellipses located at the active Gaussian centers

${(qj,pj)}j∈A$⁠. The outer rectangular boundary of each panel indicates the phase space area corresponding to the Hilbert space

$\mathcal {H}$

$H$ spanned, respectively, by the Fourier spectral and pseudospectral bases as well as the full pvN and the full pvb bases. The wavepacket, initially concentrated at the atomic core [Fig. 2(a)], is ionized by the NIR and XUV laser pulses and spreads into parts of the first (x > 0 and p > 0) and third (x < 0 and p < 0) quadrants [Figs. 2(b) and 2(c)], but a large area is never accessed. This can be intuitively expected from the classical mechanics, and indeed we see that the wavepacket closely follows the so-called simple-man trajectories⁵ (dots and + marks in Fig. 2), which obey the classical Hamiltonian of the same form as Eq. (9) with H₀ replaced by p²/2μ. In fact, we chose

$\mathcal {A}$

$A$ so that the active phase space domain contains these simple-man trajectories (with an additional margin).

In Fig. 3, we compare the photoelectron momentum distributions obtained using the reduced basis of Fig. 2 and the full pvb basis. The excellent agreement between the two results indicates that the present method not only preserves the qualitative features – cut-offs of the direct and rescattered NIR photoelectrons, and the NIR-streaked single-XUV-photon ionization peaks – but also has quantitative accuracy.

The accuracy and the size of the active pvb basis set are expected to be inversely related. However, it is not straightforward to determine the minimal size of the active set that maintains a given accuracy. To seek an upper bound to such an optimal basis size, we carried out simulations with different active sets generated by changing the margins around the simple-man trajectories. The error of a simulation was measured by ε ≔ |⟨Ψ^reduc(t_max)|Ψ^full(t_max)⟩ − 1|, where |Ψ^reduc(t_max)⟩ and |Ψ^full(t_max)⟩ are the final states calculated using a reduced and the full pvb bases, respectively. Figure 4 shows the dependence of this error on the basis set size. As the size

$NA$ is time-dependent, and as the computational cost of wavepacket propagation depends on

$N_{\mathcal {A}}$

$NA$ quadratically, we characterized

$N_{\mathcal {A}}$

$NA$ by its root-mean-square,

$\langle N_{\mathcal {A}} \rangle := \sqrt{\int _{t_{\text{min}}}^{t_{\text{max}}} dt [N_{\mathcal {A}}(t)]^2/(t_{\text{max}}-t_{\text{min}})}$

$⟨NA⟩:=∫tmintmaxdt[NA(t)]2/(tmax−tmin)$⁠, as well as its minimum and maximum values. It can be seen that the pvb basis can be compressed down to

$\langle N_{\mathcal {A}} \rangle /N = 0.14$

$⟨NA⟩/N=0.14$ or possibly further while maintaining the accuracy level at ε = 4 × 10⁻¹⁰, and at least to

$\langle N_{\mathcal {A}} \rangle /N = 0.08$

$⟨NA⟩/N=0.08$ for ε = 6 × 10⁻⁸.

The computational cost of the present method is dominated by the multiplication of the matrix

$Gjm:=(−i/ℏ)∑l∈A(Ω−1)jl⟨b̃l|H(t)|b̃m⟩$ and the vector c_m (⁠

$j,m\in \mathcal {A}$

$j,m∈A$⁠) which scales as

$O(N_{\mathcal {A}}^2)$

$O(NA2)$⁠. There is also initial overhead that scales as O(N³) originating from the computation of S_mj and (S⁻¹)_jl (m, j = 1, …, N;

$l \in \mathcal {A}$

$l∈A$⁠).

Our method can be extended straightforwardly to systems with multiple DOF.²¹ Defining the number d of DOF, the total number M of FG points increases exponentially as M = O(N^d). In contrast, the number

$MA$ of active pvb basis states for atomic and molecular systems in a laser field is expected to scale more slowly than M since the phase space coordinates are generally correlated for such systems.^4,29,30 The cost per time step of the present method goes as

$O(M_{\mathcal {A}}^2)$

$O(MA2)$⁠. The initial overhead in calculating S_mj and (S⁻¹)_jl stays at O(N³) because the multi-D Gaussians can be factored into 1D Gaussians.²¹ The dominant part in the overhead is now the calculation of the matrix elements for the potential energy operator in the reduced pvb basis, and this goes as

$O(M M_{\mathcal {A}}^2)$

$O(MMA2)$⁠. For sufficiently large d, the present method has the potential to perform better than the popular alternate-direction Crank-Nicolson scheme^31,32 [which costs O(M) per step] and the FG split-operator method^24,25 [which costs O(Mlog₂N) per step]. The pvb basis may also be used for the one-body orbitals in MCTDHF as an alternative to the static scaling of the position coordinates^8,33,34 used previously.¹⁰ 

In summary, we presented a new method to solve the TDSE based on the pvb basis. Although the basis is time-independent, the active subset is chosen in a time-dependent manner. As a first demonstration, we calculated the electronic wavepacket of a 1D atom in the combined fields of intense NIR and attosecond XUV laser pulses. This example demonstrates the high accuracy and efficiency of the method. We are currently working to extend the method to 3D with the aim of ultimately applying it to multi-electron systems in intense and ultrashort laser pulses.

This work was supported by the Minerva Foundation and the Israel Science Foundation (ISF) under Grant No. . This research is made possible by the historic generosity of the Harold Perlman family.