Long-time solution of the time-dependent Schrödinger equation for an atom in an electromagnetic field using complex coordinate contours

We demonstrate that exterior complex scaling (ECS) can be used to impose outgoing wave boundary conditions exactly on solutions of the time-dependent Schrödinger equation for atoms in intense electromagnetic pulses using finite grid methods. The procedure is formally exact when applied in the appropriate gauge and is demonstrated in a calculation of high-harmonic generation in which multiphoton resonances are seen for long pulse durations. However, we also demonstrate that while the application of ECS in this way is formally exact, numerical error can appear for long-time propagations that can only be controlled by extending the finite grid. A mathematical analysis of the origins of that numerical error, illustrated with an analytically solvable model, is also given.

Figure 1

(Color online) Upper: example exterior complex scaling contour $C\left(z\right)$ for one-dimensional problem with $z$ in bohr. Lower: propagation of a Volkov packet with mass of an electron and an initial drift momentum to the right as it oscillates across the ECS boundary at $Z_0=60\text{bohr}$. Solid lines: modulus squared of wavepacket $\left({\text{bohr}}^{-1}\right)$ on real portion of contour at three successive times, moving to the right, to the left, and again to the right. Dashed lines: wavepacket on complex portion of contour, plotted for real $z$.

Figure 2

Modulus squared of Volkov packet on ECS contour as a function of position, $z$, in bohr and time in atomic time units $\left(24.19\times {10}^{-18}\text{s}\right)$. Left column: ECS contour with $Z_0=15$, Right column: $Z_0=40$. Top row: results from method of He et al. (Ref. [21]). Middle row: correct ECS propagation. Bottom row: analytic solution of time-dependent Schrödinger equation. Quiver radius is 35.1 bohr.

Figure 3

${\left|\Psi \left(C\left(z\right),t\right)\right|}^2$ as a function of time for soft coulomb potential case. Units as in Fig. 2. Top left: results from method of He et al. (Ref. [21]) with $Z_0=25\text{bohr}$. Top right: correct ECS propagation with $Z_0=25\text{bohr}$. Quiver radius is 30.8 bohr. Bottom left: correct ECS propagation with $Z_0=60\text{bohr}$. Bottom right: propagation with real grid extending to $\pm 300\text{bohr}$ and no complex scaling.

Figure 4

Long-time propagations. Units as in Fig. 2. Left column: modulus squared of wavepackets from the soft Coulomb potential problem on ECS contour with $Z_0=60\text{bohr}$. Right column: modulus squared of Volkov state wavepackets with $Z_0=40\text{bohr}$. Top row: complex portion of ECS contour is 25 bohr in length. Middle row: complex portion is 50 bohr. Bottom row: complex portion is 75 bohr.

Figure 5

Spectrum of the time-independent operator of the alternating dc model. Operator is discretized with finite differences on the ECS contour using 100 interior points and 50 points on each of the exterior scaled domains. Spectrum of the operator is shown with circles (○); each mode is stable. The boundary of the numerical range (dashed line) significantly protrudes into the region with positive real parts, resulting in short-term growth. Dotted line: the boundary of the numerical range generated by the complex part of the ECS contour only. The rightmost point, indicated with the arrow, is the numerical abscissa of this subset of $\mathbb{C}$ and gives an estimate for the error growth discussed in Sec. 4E.

Figure 6

Time evolution of the norm of the exponential operator for the model problem with $a_0=0.6$ discretized on a finite difference grid and an ECS angle of $30°$, with time in atomic units $\left(24.19\times {10}^{-18}\text{s}\right)$. Since the operator is non-normal the numerical and spectral abscissa differ. The spectral abscissa is negative while the numerical abscissa is positive, leading to a short-term growth, but there is long-term decay. Dashed line shows the initial growth while the dotted line is the asymptotic decay rate given by the eigenvalue with the largest real part. Note that at $t=0$, the norm of $\Vert \text{exp}\left(-iHt\right)\Vert =1$.

Figure 7

For the eigenmodes for the model problem with constant vector potential of Sec. 4D1 with $\left|{\lambda }_j\right|⪡a_0$, the wave functions, $\Psi$, as functions of $z$ (in bohr) for both the fundamental modes (broken curves), are exponentially increasing on the left complex contour. The Dirichlet boundary condition forces the sum of the two to be zero at $-Z_1$. The sum of the two, shown as the solid line, is increasing in the initial part of the complex contour and then bends over to fit the homogeneous Dirichlet boundary condition at $z=-Z_1$.

Figure 8

Time evolution of total probability on interior of the domain for model problem with alternating dc field with time in atomic units $\left(24.19\times {10}^{-18}\text{s}\right)$. Solid line: exact solution based on Eq. (26); dashed line: numerical solution; dashed-dot line: error in the total probability. At time $t=400$, i.e., five periods, the error exceeds the numerical solution. Note that the error is mainly low frequency. The period is $T=80$, $a_0=0.6$, $L=\left[-13,13\right]$; we used 100 grid points on the interior of the domain. Each complex contour has a length 13 and has 50 grid points. Crank-Nicholson is used for the numerical time integration.

Figure 9

(Color online) Left column: spectral energy density in atomic units ($\text{hartrees}\times \text{atomic}$ time unit) for harmonic generation in the soft Coulomb potential problem, with 60 cycle pulse and intensities varying (top to bottom) $5\times {10}^{12}$, ${10}^{13}$, and ${10}^{14}\text{W}/{\text{cm}}^2$. Dashed vertical line: classical cutoff at $\left(I_p+3.17U_p\right)/{\omega }_0$. Right column: corresponding dipole acceleration, with time in atomic units $\left(24.19\times {10}^{-18}\text{s}\right)$.

Figure 10

Harmonic generation. Left: $S\left(\omega \right)$ from ECS calculation using a 60 cycle pulse with central frequency ${\omega }_0=0.057$ atomic units and intensity of ${10}^{14}\text{W}/{\text{cm}}^2$. Right: calculation using same pulse parameters using method of Ref. [21]. Units as in Fig. 9