Explicit schemes for time propagating many-body wave functions

Accurate theoretical data on many time-dependent processes in atomic and molecular physics and in chemistry require the direct numerical ab initio solution of the time-dependent Schrödinger equation, thereby motivating the development of very efficient time propagators. These usually involve the solution of very large systems of first-order differential equations that are characterized by a high degree of stiffness. In this contribution, we analyze and compare the performance of the explicit one-step algorithms of Fatunla and Arnoldi. Both algorithms have exactly the same stability function, therefore sharing the same stability properties that turn out to be optimum. Their respective accuracy, however, differs significantly and depends on the physical situation involved. In order to test this accuracy, we use a predictor-corrector scheme in which the predictor is either Fatunla's or Arnoldi's algorithm and the corrector, a fully implicit four-stage Radau IIA method of order 7. In this contribution, we consider two physical processes. The first one is the ionization of an atomic system by a short and intense electromagnetic pulse; the atomic systems include a one-dimensional Gaussian model potential as well as atomic hydrogen and helium, both in full dimensionality. The second process is the decoherence of two-electron quantum states when a time-independent perturbation is applied to a planar two-electron quantum dot where both electrons are confined in an anharmonic potential. Even though the Hamiltonian of this system is time independent the corresponding differential equation shows a striking stiffness which makes the time integration extremely difficult. In the case of the one-dimensional Gaussian potential we discuss in detail the possibility of monitoring the time step for both explicit algorithms. In the other physical situations that are much more demanding in term of computations, we show that the accuracy of both algorithms depends strongly on the degree of stiffness of the problem.

Figure 1

Evolution of the time step in Fatunla's method (blue line) for our Gaussian model problem. The cosine square pulse envelope (red line) is also shown on an arbitrary scale. The Gaussian potential parameters are $V_0=1$ a.u. and $\beta =1$ a.u. and we use an electromagnetic pulse with $I={10}^{14}$ W/cm${}^2$ peak intensity, $\omega =0.7$ a.u. photon energy, and a duration of 10 optical cycles.

Figure 2

Absolute error in the norm $\Delta =\left|\right(\parallel \Psi (\mathbf{r},t)\parallel -\parallel \Psi (\mathbf{r},t_0)\parallel \left)\right|$ on a logarithmic scale for different lower and upper bounds of the truncation error. The parameters of the Gaussian model problem are the same as in Fig. 1. The inset is a blow-up of the region at the beginning of the time propagation.

Figure 3

Number of Krylov vectors required to obtain convergence of the final vector propagated for different values of the (fixed) time step. The parameters of the Gaussian model problem are the same as in Fig. 1.

Figure 4

Number of iterations in the BCGSTAB and their multiplicities during the propagation with the electromagnetic pulse. The parameters of the Gaussian model problem are $V_0=4$ a.u., $\beta =0.1$ a.u., with a pulse of peak intensity $I={10}^{16}$ W/cm${}^2$, photon energy $\omega =0.5$ a.u., and duration of eight optical cycles.

Figure 5

Time-step evolution for PC method of time propagation during the propagation with the electromagnetic pulse. The parameters of the Gaussian model problem are the same as in Fig. 4.

Figure 6

Energy distribution for the Gaussian model potential. The time propagation uses Fatunla's propagator with adaptive time step and Arnoldi's propagator with five Krylov vectors and a fixed time step $\delta t=0.3$ a.u. The parameters of the model problem are $V_0=1$ and $\beta =1$ with a pulse of peak intensity I = ${10}^{14}$ W/cm${}^2$, photon energy $\omega =0.7$ a.u., and a duration of 10 optical cycles. The relative difference between both curves is of the order of ${10}^{-3}$.

Figure 7

Energy distribution for the Gaussian model potential. The time propagation uses Fatunla's propagator with adaptive time step and Arnoldi's propagator with 20 Krylov vectors and a fixed time step $\delta t=0.1$ a.u. The parameters of the model problem are $V_0=1$ and $\beta =1$ with a pulse of peak intensity $I={10}^{14}$ W/cm${}^2$, photon energy $\omega =0.1$ a.u., and duration of 10 optical cycles. The relative difference between both curves is of the order of ${10}^{-3}$.

Figure 8

Energy distribution for the Gaussian model potential. The time propagation uses Fatunla's propagator with adaptive time step and the predictor-corrector method. The parameters of the model problem are $V_0=4.0$ and $\beta =0.1$ with a pulse of peak intensity $I={10}^{16}$ W/cm${}^2$, photon energy $\omega =0.5$ a.u., and duration of eight optical cycles.

Figure 9

Energy distribution for the Gaussian model potential. The time propagation uses Arnoldi's propagator with fixed time step and the PC method with adaptive time step. The parameters of the Gaussian model problem are as in Fig. 8.

Figure 10

Energy distribution resulting from the interaction of the hydrogen atom with a cosine square pulse. Fatunla's and Arnoldi's propagators are used. The pulse has a peak intensity $I={10}^{14}$ W/cm${}^2$, a frequency $\omega =0.7$ a.u., and a duration of 10 optical cycles. The basis set of functions used is a set of 100 Coulomb Sturmian functions per angular momentum. Ten angular momenta are included and the nonlinear parameter $\kappa$ of the Coulomb Sturmian functions is equal to 0.3. The Arnoldi propagator uses five Krylov vectors and a time step of $\delta t=0.05$ a.u. The relative difference between both curves is of the order of ${10}^{-3}$.

Figure 11

Energy distribution resulting from the interaction of the hydrogen atom with a cosine square pulse. Arnoldi's propagator is used. The pulse parameters are as in Fig. 10, using 100 Coulomb Sturmian functions per angular momentum. Ten angular momenta are included and the nonlinear parameter $\kappa$ of the Coulomb Sturmian functions is equal to 0.3. For a time step of $\delta t=0.05$ a.u., we compare results when five and four Krylov vectors are used.

Figure 12

Energy distribution resulting from the interaction of the hydrogen atom with a cosine square pulse. Arnoldi's propagator is used. The pulse has a peak intensity I = 10${}^{14}$ W/cm${}^2$, a frequency $\omega =0.114$ a.u., and a duration of 20 optical cycles. We use a set of 600 Coulomb Sturmian functions per angular momentum. Ten angular momenta are taken into account and the nonlinear parameter of the Coulomb Sturmian functions $\kappa =0.3$. The Arnoldi propagator uses 25 Krylov vectors and a time step of $\delta t=0.05$ a.u. The relative difference between both curves is of the order of ${10}^{-3}$.

Figure 13

Energy distribution resulting from the interaction of the hydrogen atom with a cosine square pulse. Arnoldi's propagator is used. The pulse has a peak intensity I = 10${}^{14}$ W/cm${}^2$, a frequency $\omega =0.0228$ a.u., and a duration of four optical cycles. The basis set of functions used is a set of 1200 Coulomb Sturmian functions per angular momentum. Eighty angular momenta are taken into account and the nonlinear parameter $\kappa$ of the Coulomb Sturmian functions is equal to 0.3. The Arnoldi propagator uses 70 Krylov vectors and a time step of $\delta t=0.05$ a.u.

Figure 14

Single ionization spectrum resulting from the interaction of a helium atom with a short cosine square pulse. Arnoldi's propagator is used. The pulse has a peak intensity of I = 10${}^{14}$ W/cm${}^2$, a frequency $\omega =2.1$ a.u., and a duration of six optical cycles. The basis set of functions used is a set of 140 B-spline functions of order seven per electron angular momentum. The total angular momentum L = 0,1,2 and the maximum value of the individual electron angular momentum is three. The box size is 150 a.u. The Arnoldi propagator uses 40 Krylov vectors.

Figure 15

Susceptibility $\chi \left(t\right)$ calculated for a quantum dot with $\omega =$1.0 a.u., using the Arnoldi's propagator. We took five Krylov vectors and compared results for two different values of the time step.

Figure 16

Susceptibility $\chi \left(t\right)$ calculated for a quantum dot with $\omega =$1.0 a.u., using the Arnoldi's propagator. We take two different combinations of time step and of the dimension of the Krylov subspace used, to illustrate the compromise between these two quantities.