Adaptive Exponential Integrators for MCTDHF

TL;DR

This paper evaluates various exponential integrators for the MCTDHF method, finding that exponential Lawson multistep methods with predictor/corrector steps offer optimal stability, accuracy, and efficiency, especially with adaptive time-stepping.

Contribution

It introduces and assesses exponential Lawson multistep integrators with predictor/corrector steps for MCTDHF, demonstrating their superior stability and efficiency over traditional methods.

Findings

Exponential Lawson multistep methods outperform other integrators in stability and accuracy.

Predictor step enables reliable adaptive time-stepping without extra cost.

Optimal balance of computational cost and accuracy achieved with these methods.

Abstract

We compare exponential-type integrators for the numerical time-propagation of the equations of motion arising in the multi-configuration time-dependent Hartree-Fock method for the approximation of the high-dimensional multi-particle Schr{\"o}dinger equation. We find that among the most widely used integrators like Runge-Kutta, exponential splitting, exponential Runge-Kutta, exponential multistep and Lawson methods, exponential Lawson multistep methods with one predictor/corrector step provide optimal stability and accuracy at the least computational cost, taking into account that the evaluation of the nonlocal potential terms is by far the computationally most expensive part of such a calculation. Moreover, the predictor step provides an estimator for the time-stepping error at no additional cost, which enables adaptive time-stepping to reliably control the accuracy of a computation.

Equations

$$\mathrm{i}{\partial }_{t}u(t)=Au(t)+B(u(t))=H(u(t)),t>t_0,u(t_0)=u_0,$$

$$\mathrm{i}\frac{\partial \psi }{\partial t}=H\psi ,$$

$$H$$

$$T$$

$$u=\sum _{(j_1,\dots ,j_f)}{a}_{j_1,\dots ,j_f}(t){\varphi }_{j_1}(x_1,t)\cdots {\varphi }_{j_f}(x_f,t)=:\sum _J{a}_J(t){\Phi }_J(x,t).$$

$$\displaystyle{\mathrm{i}}\,\frac{da_{J}}{dt}=\sum_{K}\big{\langle}\Phi_{J}\left|W\right|\Phi_{K}\big{\rangle}\,a_{K}\quad\forall J,$$

$$\mathrm{i}\frac{\partial {\varphi }_j}{\partial t}=T{\varphi }_j+(I-P)\sum _{k=1}^N\sum _{\ell =1}^N{\rho }_{j,\ell }^{-1}{\overline{W}}_{\ell ,k}{\varphi }_k,j=1,\dots ,N,$$

$$\overline{W}_{j,\ell}=\big{\langle}\psi_{j}|W|\,\psi_{\ell}\big{\rangle},\quad\mbox{with}\quad\psi_{j}=\big{\langle}\phi_{j}\,|\,u\big{\rangle},~{}~{}\rho_{j,\ell}=\big{\langle}\psi_{j}\,|\,\psi_{\ell}\big{\rangle},$$

$$A=\frac{{\mathrm{i}}}{2}\Big{(}0,\Delta^{(1)},\dots,0,\Delta^{(f)}\Big{)}^{T},\quad B=B(a,\phi),$$

$$u_{n+1}:=\mathcal{S}(h,u_n)={\mathcal{E}}_B(b_{s}h,\cdot )\circ {\mathcal{E}}_A(a_{s}h,\cdot )\circ \dots \circ {\mathcal{E}}_B(b_{1}h,\cdot )\circ {\mathcal{E}}_A(a_{1}h,u_n),$$

$$u(t_n+h)={\mathrm{e}}^{hA}{u}_n+{\int }_0^h{\mathrm{e}}^{(h-\tau )A}B(u(t_n+\tau ))\mathrm{d}\tau .$$

$$\big{(}t_{n}+h\tau_{1},B(u(t_{n}+h\tau_{1})\big{)},\dots,\big{(}(t_{n}+h\tau_{k},B(u(t_{n}+h\tau_{k})\big{)}.$$

$$\big{(}-(k-1)\,h,B(u_{n-k+1})\big{)},\dots,\big{(}(-h,B(u_{n-1})),(0,B(u_{n})\big{)}\,.$$

$$u^{\prime }(t)={\mathrm{e}}^{-tA}B({\mathrm{e}}^{tA}u(t))=:F(u(t)),$$

$$H(t)=H_0+(x_1+x_2)\mathcal{E}(t),$$

$$H_0=-\frac{1}{2}({\partial }_{x_1}^2+{\partial }_{x_2}^2)-\frac{2}{\sqrt{x_1^2+b^2}}-\frac{2}{\sqrt{x_2^2+b^2}}+\frac{1}{\sqrt{(x_1-x_2{)}^2+b^2}},$$

$$\mathcal{E}(t)={\mathcal{E}}_{0}g(t)\sin (\omega t).$$

Full text

1 Introduction

We compare time integration methods for nonlinear Schrödinger-type equations

[TABLE]

on the Hilbert space $\mathcal{B}{{=L^{2}}}$. Here, $A\colon\mathcal{D}\subseteq\mathcal{B}\to\mathcal{B}$ is a self-adjoint differential operator and $B$ a generally unbounded nonlinear operator. Our focus is on the equations of motion associated with the multi-configuration time-dependent Hartree-Fock (MCTDHF) approximation to the multi-particle electronic Schrödinger equation, where the key issue is the high computational effort for the evaluation of the nonlocal (integral) operator $B$. Thus, in the choice of the most appropriate integrator, we emphasize a minimal number of evaluations of $B$ for a given order and disregard the effort for the propagation of $A$, which can commonly be realized at essentially the cost of two (cheap) transforms between real and frequency space via fast transforms like [I]FFT. The approaches that we pursue and advocate in this paper are thus based on splitting of the vector fields in (1). It turns out that exponential integrators hocost10 based on the variation of constants serve our purpose best, as they provide a desirable balance between computational effort and stability.

2 The MCTDHF method

We focus on the comparison of numerical methods for the equations of motion associated with MCTDHF for the approximate solution of the time-dependent multi-particle Schrödinger equation

[TABLE]

where the complex-valued wave function $\psi=\psi(x_{1},\dots,x_{f},t)$ depends on time $t$ and, in the case considered here, the positions $x_{1},\dots,x_{f}\in\mathbf{R}^{3}$ of electrons in an atom or molecule. The time-dependent Hamiltonian reads

[TABLE]

Here $V_{\mathrm{ext}}(x_{1},\dots,x_{f},t)$ is a smooth time-dependent function, and $\Delta^{(k)}$ is the Laplace operator with respect to $x_{k}$ only.

In MCTDHF as put forward in zanghellinietal03c , the multi-electron wave function $\psi$ is approximated by a function $u$ living in a manifold $\mathcal{M}$ characterized by the ansatz

[TABLE]

For the electronic Schrödinger equation, the Pauli principle implies that only solutions $u$ are considered which are antisymmetric under exchange of any pair of arguments $x_{j},x_{k}$,

Now, the Dirac-Frenkel variational principle frenkel34 in conjunction with orthogonality conditions is used to derive differential equations for the coefficients $a_{J}$ and the so-called single-particle functions $\phi_{j}$ in (2), where we will henceforth tacitly identify $u$ with the vector $(a,\phi)$ of coefficients and orbitals,

[TABLE]

where

[TABLE]

and where $P$ is the orthogonal projector onto the space spanned by the functions $\phi_{j}$. We will henceforth denote

[TABLE]

where $B$ is the vector of the components associated with the potential which constitute the computationally most expensive part.

2.1 Splitting methods

Popular integrators for quantum dynamics are exponential time-splitting methods which are based on multiplicative combinations of the partial flows $\mathcal{E}_{A}(t,u)\colon$ $u\mapsto u(t)={\mathrm{e}}^{tA}\,u$ and $\mathcal{E}_{B}(t,u)\colon u\mapsto u(t)$ with $u^{\prime}(t)=B(u(t)),\;u(0)=u$. For a single step $(t_{n},u_{n})\mapsto(t_{n+1},u_{n+1})$ with time-step $h$, this reads

[TABLE]

where the coefficients $a_{j},b_{j},\,j=1\ldots s$ are determined according to the requirement that a prescribed order of consistency is obtained haireretal02 . For a convergence analysis of splitting methods in the context of MCTDHF, see for instance knth10a .

2.2 Exponential integrators

An approach which also exploits the separated vector fields is given by the class of exponential integrators, which are comprehensively discussed in hocost10 . Here the variation of constant formula is used to express the solution of (1) for a time-step $t_{n}\to t_{n+1}=t_{n}+h$ via the integral equation

[TABLE]

Different numerical integrators are distinguished depending on how the integral in (6) is approximated.

Exponential Runge-Kutta methods

When the integral in (6) is approximated by a quadrature formula of Runge-Kutta type, relying on evaluations of the nonlinear operator $B$ at interior points $t_{n}+h\tau_{j},\ \tau_{j}\in[0,1],\ j=1,\dots,k$, an exponential Runge-Kutta method is obtained. This corresponds to replacing $B(\cdot)$ in the integrand by a polynomial interpolant at the points

[TABLE]

The method is realized by stepping from $t_{n}+h\tau_{j}\to t_{n}+h\tau_{j+1}$ in the same way as for a Runge-Kutta method, with appropriate weights of the stages. For implicit methods, nonlinear systems of equations have to be solved, which is generally considered as prohibitive. Note that after interpolation, the resulting integral can be evaluated analytically by using the $\varphi$-functions or alternatively, by numerical quadrature hocost10 . Such a procedure has first been proposed in friedli78 , for a stiff error analysis, see hocost10 and references therein. For our comparisons, we use the fourth order Krogstad method mentioned there.

Exponential multistep methods

The integral in (6) can be approximated in terms of an interpolation polynomial at previous approximations

[TABLE]

This yields an (explicit) exponential Adams-Bashforth multistep method first mentioned in certaine60 , and introduced more systematically in noersett69 , see also for instance calvopal11 ; hocost10 . If the interpolation also comprises the forward point $(h,B(u_{n+1}))$, an (implicit) exponential Adams-Moulton method is obtained. These two approaches can be combined in a predictor/corrector method in the same way as for linear multistep methods. Exponential multistep methods have first been considered and analyzed in calvopal11 under the assumption of smooth $B$, and a starting strategy is also given there.

Lawson methods

In Lawson methods, equation (1) is transformed prior to the numerical integration by the substitution $u(t)\to{\mathrm{e}}^{-tA}\,u(t)$. To the resulting equation

[TABLE]

any appropriate time-stepping scheme can be applied. The main advantage lies in the fact that the dynamics associated with the non-smooth operator $A$ is separated by the transformation which can be realized cheaply in frequency space, while the problem subjected to the time-stepping scheme is smoother, thus allowing for larger time-steps. This transformation was first introduced in lawson67 for ordinary differential equations.

In a one-step version, an explicit Runge–Kutta method is employed to solve (8), which is equivalent to interpolation at interior nodes of the whole integrand in (6) by a polynomial in the same fashion as in (7). Reference hocost17 gives a convergence proof of Lawson Runge-Kutta methods in the stiff case, however under the assumption that the operator $B$ is smooth, which is not the case in the MCTDHF equations we are considering. A convergence proof for Adams-Lawson multistep methods for the MCTDHF equations under minimal regularity requirements is given in the forthcoming work koch19 . The proof addresses the transformed equation (8) and combines stability and consistency to conclude convergence. To this end, a boot-strapping argument is employed, first showing convergence in the Sobolev space $H^{1}$. Stability in $L^{2}$ only holds if the numerical solution is in $H^{1}$, which follows from the first argument, whence convergence in $L^{2}$ is inferred. Lipschitz conditions for the right-hand side entering the stability arguments can be shown by appropriate Sobolev-type inequalities in both $H^{1}$ and $L^{2}$. To prove consistency, the norms of derivatives of $F$ in (8) are estimated, which amounts to bounds on commutators of the operators $A$ and $B$. This implies assumptions on the regularity of the exact solution $u$.

We will demonstrate that the best approach for our goal is to use exponential Lawson multistep methods in a predictor/corrector implementation, which is shown to increase the accuracy and also provides a local error estimator for adaptive time-stepping at no additional cost. The efficiency of the time discretization can be improved if high-order time propagators are employed. In the multistep approach, this does not imply additional computational cost if no memory limitations have to be taken into account.

Comparisons

To assess the performance of the exponential integration methods described above, we will also show results for the classical explicit Runge-Kutta method of fourth order (RK4) and the second-order Strang splitting.

3 Numerical results

To illustrate the performance of our numerical methods, we consider MCTDHF with the choice $N=4$ for a one-dimensional model of a helium atom investigated in zanghellinietal03c , where111Note that in exponential integrators, the explicit time-dependence in the potential does not call for a special treatment in the numerical quadrature, in the splitting methods, the potential is propagated by an explicit Runge-Kutta method of appropriate order.

[TABLE]

with a smoothed Coulomb potential with shielding parameter $b=0.7408$, which is irradiated by a short, intense, linearly polarized laser pulse

[TABLE]

The peak amplitude is set to $\mathcal{E}_{0}=0.1894$, the frequency is $\omega=0.1837$, and we define the envelope g(t)=1.2\exp\big{(}-5\cdot 10^{-4}\left(t-6\pi/\omega\right)^{2}\big{)}. The parameters are taken from zanghellinietal03c , and the envelope is a smooth approximation of the trapezoidal envelope chosen there. In zanghellinietal03c , this model serves to illustrate the effect of correlation on the probability density along the diagonal $x=y$, which implies that the single-configuration Hartree–Fock approximation is insufficient. We first investigate stable long-time propagation in Fig. 1. We monitor norm conservation of the wave function in the propagation of the ground state for the Hamiltonian $H_{0}=H(0)$ for different equidistant stepsizes to resolve precisely the onset of instability. For RK4, the number of steps is specified in the plot; for all other methods, the number of steps is in $\{1000,2000,\dots,12000\}$. If norm conservation is violated beyond the effect of numerical accuracy, the method cannot be recommended for physical applications. Indeed, we observe the following: Explicit Runge-Kutta methods only behave in a stable way when the numerical accuracy is already very high, close to round-off error. Exponential multistep methods222In this experiment, all multistep methods are started by the Krogstad exponential Runge-Kutta method with stepsize $h/50$. behave stably only for short times, even when a corrector step is performed. Exponential Runge-Kutta and Runge-Kutta-Lawson methods behave stably, likewise as splitting methods. Adams-Lawson multistep methods behave very stably, a corrector step adds to the accuracy, as well as providing an error estimate as the basis for adaptive time-stepping. The unstable exponential multistep methods are no longer considered. While showing the same stability behavior, the Yoshida splitting is demonstrated to be less efficient than the Suzuki splitting, and the low order (but popular) Strang splitting is not competitive. Higher-order multistep methods provide higher accuracy at the same computational effort irrespective of the order, and are thus also considered for this comparison.

Next, we compare the efficiency of the different integrators. The unstable exponential multistep methods are no longer considered. While showing the same stability behavior, the Yoshida splitting is demonstrated to be less efficient than the Suzuki splitting, and the low order (but popular) Strang splitting is not competitive. High-order multistep methods provide higher accuracy at the same computational effort and are thus also considered for this comparison. To this end, we plot in Fig. 2 the accuracy as compared to a very precise reference solution at $t=80$ as a function of the number of evaluations of the computationally expensive potential part $B$ (dots on solid lines). Furthermore, we give the CPU time required in a sequential implementation on one thread of the Vienna Scientific Cluster (VSC) 3 comprising one Intel Xeon E5-2650v2 processor with 8 kernels of 2.6 gHz (crosses ‘$\times$’). We note that, as expected, the runtime is proportional to the number of potential evaluations. We observe that high-order Lawson multistep methods perform best, where particularly the high order which can be achieved in the multistep versions without additional evaluations is advantageous. Splitting methods, particularly the low order Strang splitting, are not very efficient due to the high effort for the propagation of the potential.

We stress that this shows only the picture on uniform grids. The multistep versions show their biggest advantage in adaptive time-stepping due to the cheap means of error estimation in the predictor/corrector implementation. To demonstrate that this works reliably for Adams-Lawson methods, we show in Fig. 3 the laser field $\mathcal{E}(t)$ and total energy functional (top) illustrating the local solution smoothness, and the stepsizes (bottom) automatically generated for the Adams-Lawson method of order 6. We see that the adaptively chosen stepsizes reflect the smoothness of the time evolution and the Lawson method enables larger stepsizes. The Adams-Lawson solution has been confirmed to be converged to within the prescribed tolerance $10^{-5}$. On the other hand, the corresponding exponential multistep method (without the Lawson transformation) shows a noticeable deviation in the solution.