Efficient solver for time-dependent Schrödinger equation with interaction between atoms and strong laser field^*

Theoretical simulation plays an important role in understanding of the interaction of atoms or molecules with a strong laser field. There are various methods of simulating the strong-field physical processes, such as the simple-man model,^[1] strong-field approximation (SFA),^[2,3] classical-trajectory Monte Carlo simulation (CTMC),^[4] quantum-trajectory Monte Carlo (QTMC),^[5,6] etc. While the most accurate theoretical method of obtaining the complete information is to solve the time-dependent Schrödinger equation (TDSE). However, solving the three-dimensional (3D) TDSE for simulating such laser interactions with atoms becomes a challenge because, with the development of laser technology, the laser field intensity can be much higher than the over-barrier ionization region of atoms or molecules. Although many efforts have been made in the past few decades^[7–14] to tackle this challenging problem, it still needs more efficient numerical methods, parallel to adapting to the development of computer technology.

Several efficient parallel methods of solving the 3D TDSE^[15–18] have been developed, mainly based on direct TDSE solutions in Cartesian coordinates. Considering the atomic spherical symmetry, it should be better to choose spherical coordinates in computation but this is difficult because, in general, the methods with using the basis set expansion techniques in spherical coordinates are not suitable for parallel calculation.^[15] To solve this problem, Patchkovskii et al. have proposed a parallel method,^[18] and achieved the spatial derivatives of the radial coordinate by using the finite-difference method (FDM), which leads to an increase in the grids for the larger sized problem,^[17] and finally significantly increasing the amount of computation in the case of higher accuracy.^[19] Guan et al. solved TDSE of two-electron quantum system in spherical coordinates efficiently based on the finite-element (FE) discrete-variable-representation (DVR) and Short-Iteration-Lanczos (SIL) propagator,^[20] because the FE-DVR can offer higher computing accuracy^[21] and make the Hamiltonian matrix very sparse, whose effects are similar to those provided by the B-spline^[22,23] and sine-DVR,^[24] and the efficiency of the wave function evolution strongly relies on the sparsity of the Hamiltonian matrix for the SIL propagator.^[25] The most ideal structure of the Hamiltonian matrix is expected to have diagonal or similarly sparse elements for the SIL propagator.

In this work, we develop a parallel numerical method of solving TDSE in spherical coordinates for simulating linearly polarized laser interacting with single active electron (SAE) atoms. First, a sparse Hamiltonian matrix is constructed based on the FE-DVR and spherical harmonics functions. Second, ideal sub-matrices are separated out from the Hamiltonian matrix for the SIL propagator,^[26,27] which is based on the real space product formula (RSPF). This combination guarantees efficient parallel computing for the simulation of real physical problems. We then apply this method to the efficient calculation of atomic high harmonic generation (HHG). Atomic units are used throughout this paper unless otherwise indicated.

To simulate the interaction of atoms with linear polarized laser fields,^[24,28] the TDSE for the full wave function can be expressed as follows:

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{\partial }{\partial t}\psi ({\boldsymbol{r}},t)=\left(-\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}+V(r)+{V}_{I}({\boldsymbol{r}},t)\right)\psi ({\boldsymbol{r}},t),\end{eqnarray} \tag{ 1 }$$

with the potential term $V(r)=-1/r$ (atom hydrogen) and the laser–matter interaction term defined as

$$\begin{eqnarray}&&{V}_{I}({\boldsymbol{r}},t)={\boldsymbol{r}}\cdot {\boldsymbol{E}}(t).\end{eqnarray} \tag{ 2 }$$

The length gauge is adopted in this computation. Thus, the Schrödinger equation for the reduced wave functions can be rewritten as follows:

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{\partial }{\partial {\rm{t}}}\psi ({\boldsymbol{r}},{\rm{t}})=\left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {{\rm{r}}}^{2}}+\displaystyle \frac{{{\rm{L}}}^{2}}{2{{\rm{r}}}^{2}}-\displaystyle \frac{1}{{\rm{r}}}+\cos \theta \mathrm{rE}({\rm{t}})\right)\psi ({\boldsymbol{r}},{\rm{t}})\end{eqnarray} \tag{ 3 }$$

with

$$\begin{eqnarray}&&{L}^{2}=-\displaystyle \frac{1}{\sin \theta }\displaystyle \frac{\partial }{\partial \theta }\left(\sin \theta \displaystyle \frac{\partial }{\partial \theta }\right)-\displaystyle \frac{1}{\sin {\theta }^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\phi }^{2}}.\end{eqnarray} \tag{ 4 }$$

To solve Eq. (3), the reduced wave function $\psi ({\boldsymbol{r}},t)$ is expanded into a series of spherical harmonics ${{\rm{Y}}}_{{lm}}(\theta,\phi )$ and radial functions ${\varphi }_{{lm}}(r,t)$,

$$\begin{eqnarray}&&\psi ({\boldsymbol{r}},t)=\displaystyle \frac{1}{r}\displaystyle \underset{l=0}{\overset{{L}_{\max }}{\sum }}\displaystyle \underset{m=-l}{\overset{l}{\sum }}{\varphi }_{{lm}}(r,t){{\rm{Y}}}_{{lm}}(\theta,\phi ).\end{eqnarray} \tag{ 5 }$$

For a linear laser field, none of photons have angular momentum, so that all calculations can be done with m = 0 if the initial atomic state has a zero magnetic quantum number. Thus, equation (5) is reduced as follows:

$$\begin{eqnarray}&&\psi ({\boldsymbol{r}},t)=\displaystyle \frac{1}{r}\displaystyle \underset{l=0}{\overset{{L}_{\max }}{\sum }}{\varphi }_{l}(r,t){{\rm{Y}}}_{l0}(\theta ).\end{eqnarray} \tag{ 6 }$$

The partial waves can be further discretized as ${\varphi }_{l}(r,t)={\sum }_{k}{c}_{k}^{l}(t){\chi }_{k}(r)$, where ${\chi }_{k}(r)$ is FE-DVR base. Two vectors are set to be $X(r)=({\chi }_{1}(r),{\chi }_{2}(r),\ldots,{\chi }_{k}(r))$ and ${C}_{l}(t)=({c}_{1}^{l}(t),{c}_{2}^{l}(t),\ldots,{c}_{k}^{l}(t))$. Substituting Eq. (6) into Eq. (3), and separating the time-dependent parts from time-independent parts of the Hamiltonian matrix, the corresponding diagonal matrix element and off-diagonal matrix element are given by^[24]

$$\begin{eqnarray}&&\begin{array}{l}{h}_{l,l}=\left\langle \displaystyle \frac{{X}^{T}(r)}{r}\right|-\displaystyle \frac{{\partial }^{2}}{2\partial {r}^{2}}+\displaystyle \frac{l(l+1)}{2{r}^{2}}-\displaystyle \frac{1}{r}\left|\displaystyle \frac{X(r)}{r}\right\rangle,\\ l=0,1,2,\ldots,\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{h}_{l,l\pm 1}(t)=\left\langle \displaystyle \frac{{X}^{T}(r)}{r}\right|{g}_{l\pm 1}{rE}(t)\left|\displaystyle \frac{X(r)}{r}\right\rangle,\,\,\,l\gt 0,\end{eqnarray} \tag{ 8 }$$

where g_l is

$$\begin{eqnarray}&&{g}_{l}=\displaystyle \frac{l+1}{\sqrt{(2l+1)(2l+3)}}.\end{eqnarray} \tag{ 9 }$$

The Hamiltonian matrix is symmetrically tridiagonal, as indicated below.

$$\begin{eqnarray}H(r,t)=\left(\begin{array}{cccccc}{h}_{\mathrm{0,0}} & {h}_{\mathrm{0,1}}(t) & O & \cdots & \cdots & O\\ {h}_{\mathrm{0,1}}(t) & {h}_{\mathrm{1,1}} & {h}_{\mathrm{2,1}}(t) & O & \cdots & O\\ O & {h}_{\mathrm{1,2}}(t) & {h}_{\mathrm{2,2}} & \ddots & \ddots & \vdots \\ \vdots & O & \ddots & \ddots & {h}_{l-1,l-2}(t) & O\\ \vdots & \vdots & \ddots & {h}_{l-2,l-1}(t) & \ddots & {h}_{l,l-1}(t)\\ O & O & \cdots & O & {h}_{l,l-1}(t) & {h}_{l,l}\end{array}\right),\end{eqnarray} \tag{ 10 }$$

where O is the zero matrix. According to the RSPF, H(r,t) can be divided into two matrices A and B(t),

$$\begin{eqnarray}A=\left(\begin{array}{cccccc}{h}_{\mathrm{0,0}} & O & O & \cdots & \cdots & O\\ O & {h}_{\mathrm{1,1}} & O & O & \cdots & O\\ O & O & {h}_{\mathrm{2,2}} & \ddots & \ddots & \vdots \\ \vdots & O & \ddots & \ddots & O & O\\ \vdots & \vdots & \ddots & O & \ddots & O\\ O & O & \cdots & O & O & {h}_{l,l}\end{array}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}B(t)=\left(\begin{array}{cccccc}O & {h}_{\mathrm{1,0}}(t) & O & \cdots & \cdots & O\\ {h}_{\mathrm{0,1}}(t) & O & {h}_{\mathrm{2,1}}(t) & O & \cdots & O\\ O & {h}_{\mathrm{1,2}}(t) & O & \ddots & \ddots & \vdots \\ \vdots & O & \ddots & \ddots & {h}_{l-1,l-2}(t) & O\\ \vdots & \vdots & \ddots & {h}_{l-2,l-1}(t) & \ddots & {h}_{l,l-1}(t)\\ O & O & \cdots & O & {h}_{l-1,l}(t) & O\end{array}\right).\end{eqnarray} \tag{ 12 }$$

Here, B(t) is not a diagonal matrix and can be reshaped further. After setting

$$\begin{eqnarray}{b}_{l}(t)=\left(\begin{array}{cc}O & {h}_{l,l-1}(t)\\ {h}_{l-1,l}(t) & O\end{array}\right),\end{eqnarray} \tag{ 13 }$$

B(t) can be decomposed into two block diagonal matrices as follows:

$$\begin{eqnarray}\begin{array}{rcl}B(t) & = & {B}_{\mathrm{odd}}(t)+{B}_{\mathrm{even}}(t)\\ & = & \left(\begin{array}{cccc}{b}_{1}(t) & O & \cdots & O\\ O & {b}_{3}(t) & \ddots & \vdots \\ \vdots & \ddots & \ddots & O\\ O & \cdots & O & \ddots \end{array}\right)+\left(\begin{array}{cccc}O & O & \cdots & O\\ O & {b}_{2}(t) & \ddots & \vdots \\ \vdots & \ddots & {b}_{4}(t) & O\\ O & \cdots & O & \ddots \end{array}\right).\end{array}\end{eqnarray} \tag{ 14 }$$

Using Lie–Trotter–Suzuki product formula (LTSPF),^[29,30] the propagation scheme can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}U(H,t,t+\delta t) & = & \exp (-{\rm{i}}\delta t(A+{B}_{\mathrm{odd}}(t)+{B}_{\mathrm{even}}(t)))\\ & = & \exp \left(-\displaystyle \frac{{\rm{i}}\delta t}{2}{B}_{\mathrm{odd}}(t)\right)\exp \left(-\displaystyle \frac{{\rm{i}}\delta t}{2}{B}_{\mathrm{even}}(t)\right)\\ & & \times \exp (-{\rm{i}}\delta \mathrm{tA})\exp \left(-\displaystyle \frac{{\rm{i}}\delta {\rm{t}}}{2}{{\rm{B}}}_{\mathrm{even}}({\rm{t}})\right)\\ & & \exp \left(-\displaystyle \frac{{\rm{i}}\delta {\rm{t}}}{2}{B}_{\mathrm{odd}}(t)\right)+O(\delta {t}^{3}).\end{array}\end{eqnarray} \tag{ 15 }$$

As shown in Eqs. (11) and (14), the matrices A, ${B}_{\mathrm{even}}(t)$, and ${B}_{\mathrm{odd}}(t)$ are block-diagonal matrices. The sub-matrices in these three matrices act independently on the corresponding partial waves. Thus, the parallelization of propagation scheme can be realized programmatically which can be used very well to solve TDSE as shown in a two-electron case.^[31] From a concrete realization of the propagation $\exp (-{\rm{i}}\delta \mathrm{tA})$, for each ${h}_{l,l}$ in matrix A,

$$\begin{eqnarray}&&{C}_{l}^{{\rm{T}}}(t+\delta t)=\exp (-{\rm{i}}\delta {{th}}_{l,l}){C}_{l}^{{\rm{T}}}(t).\end{eqnarray} \tag{ 16 }$$

And for each ${b}_{l}(t)$ in B_odd and B_even,

$$\begin{eqnarray}&&\begin{array}{l}\left(\begin{array}{l}{C}_{l-1}^{{\rm{T}}}(t+\delta t)\\ {C}_{l}^{{\rm{T}}}(t+\delta t)\end{array}\right)\\ =\exp \left(\displaystyle \frac{-{\rm{i}}\delta t}{2}{b}_{l}(t)\right)\left(\begin{array}{l}{C}_{l-1}^{{\rm{T}}}(t)\\ {C}_{l}^{{\rm{T}}}(t)\end{array}\right)\\ =\exp \left(\displaystyle \frac{-{\rm{i}}\delta {\rm{t}}}{2}\left[\begin{array}{cc}O & {h}_{l,l-1}(t)\\ {h}_{l-1,l}(t) & O\end{array}\right]\right)\left(\begin{array}{l}{C}_{l-1}^{{\rm{T}}}(t)\\ {C}_{l}^{{\rm{T}}}(t)\end{array}\right).\end{array}\end{eqnarray} \tag{ 17 }$$

According to the space division scheme of FE-DVR,^[32] the radial coordinate of the interval $[0,{r}_{\max }]$ is divided into n_e finite elements. In each FE ($i.e.$, in $[{r}_{i},{r}_{i+1}])$, n_g local DVR bases ${\chi }_{i,m}$ are constructed based on the generalized Gauss–Lobatto points r_im and weights w_im. In the numerical computation, ${\chi }_{k}(r)={\chi }_{i,m}(r)$. So, the partial waves can be expressed as ${\varphi }_{l}(r,t)={\sum }_{i,m}{c}_{{im}}^{l}(t){\chi }_{i,m}(r)$. The boundary conditions that the wave function vanishes at ${r}_{11}=0$ and ${r}_{{n}_{{\rm{e}}}{n}_{{\rm{g}}}}={r}_{\max }$ are imposed by deleting the first and last functions ${\chi }_{\mathrm{1,1}}$ and ${\chi }_{{n}_{{\rm{e}}},{n}_{{\rm{g}}}}$. There are ${n}_{\mathrm{basis}}={n}_{{\rm{e}}}{n}_{{\rm{g}}}-{n}_{{\rm{g}}}-1$ basis functions in use.

According to the characteristics of the FE-DVR bases, the potential-energy matrix and the matrix of any other local operator turn out to be diagonal with regard to the elements i and local DVR basis indices m; i.e.,

$$\begin{eqnarray}&&\begin{array}{l}\left\langle \displaystyle \frac{{\chi }_{i,m}(r)}{r}\right|{V}_{l}(r,t)\left|\displaystyle \frac{{\chi }_{i^{\prime},m^{\prime} }(r)}{r}\right\rangle \\ =\left\langle \displaystyle \frac{{\chi }_{i,m}(r)}{r}\right|\displaystyle \frac{l(l+1)}{2{r}^{2}}-\displaystyle \frac{1}{r}\left|\displaystyle \frac{{\chi }_{i^{\prime},m^{\prime} }(r)}{r}\right\rangle \\ ={\delta }_{{ii}^{\prime} }{\delta }_{{mm}^{\prime} }{V}_{l}({r}_{i,m})={v}_{{im},{im}}^{l},\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}{V}_{l}=\left(\begin{array}{cccccccc}{v}_{11,11}^{l} & O & 0 & O & 0 & \cdots & \cdots & O\\ O & \ddots & O & O & O & \ddots & & \vdots \\ 0 & O & {v}_{21,21}^{l} & O & 0 & O & \ddots & O\\ O & O & O & \ddots & O & O & O & \cdots \\ 0 & O & 0 & O & {v}_{31,31}^{l} & O & 0 & \cdots \\ \vdots & \ddots & O & O & O & \ddots & O & \vdots \\ \vdots & & \vdots & O & 0 & O & {v}_{41,41}^{l} & O\\ O & \cdots & O & \vdots & \vdots & \cdots & \cdots & \ddots \end{array}\right).\end{eqnarray} \tag{ 19 }$$

While for the kinetic-energy representation, the matrix is shown in a diagonal overlapping block structure as follows:^[33]

$$\begin{eqnarray}&&{t}_{{im},{im}^{\prime} }=\left\langle \displaystyle \frac{{\chi }_{i,m}(r)}{r}\right|-\displaystyle \frac{1}{2}{{\rm{\nabla }}}^{2}\left|\displaystyle \frac{{\chi }_{i^{\prime},m^{\prime} }(r)}{r}\right\rangle r,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}T=\left(\begin{array}{cccccccc}{t}_{\mathrm{11,11}} & \cdots & {t}_{1m,11} & O & 0 & \cdots & \cdots & O\\ \vdots & \ddots & \vdots & O & O & \ddots & & \vdots \\ {t}_{\mathrm{11,1}m} & \cdots & {t}_{\mathrm{21,21}} & \cdots & {t}_{2m,21} & O & \ddots & O\\ O & O & \vdots & \ddots & \vdots & O & O & \cdots \\ 0 & O & {t}_{\mathrm{21,2}m} & \cdots & {t}_{\mathrm{31,31}} & \cdots & {t}_{3m,31} & \cdots \\ \vdots & \ddots & O & O & \vdots & \ddots & \vdots & \vdots \\ \vdots & & \vdots & O & {t}_{\mathrm{31,3}m} & \cdots & {t}_{3m,3m} & \vdots \\ O & \cdots & O & \vdots & \vdots & \cdots & \cdots & \ddots \end{array}\right).\end{eqnarray} \tag{ 21 }$$

Therefore, the elements of matrix ${h}_{l,l}$ in Eq. (11) can be expressed as

$$\begin{eqnarray}{h}_{l,l}=T+{V}_{l}=\left(\begin{array}{cccccccc}{t}_{\mathrm{11,11}}+{v}_{11,11}^{l} & \cdots & {t}_{1m,11} & O & 0 & \cdots & \cdots & O\\ \vdots & \ddots & \vdots & O & O & \ddots & & \vdots \\ {t}_{\mathrm{11,1}m} & \cdots & {t}_{\mathrm{21,21}}+{v}_{21,21}^{l} & \cdots & {t}_{2m,21} & O & \ddots & O\\ O & O & \vdots & \ddots & \vdots & O & O & \cdots \\ 0 & O & {t}_{\mathrm{21,2}m} & \cdots & {t}_{\mathrm{31,31}}+{v}_{31,31}^{l} & \cdots & {t}_{3m,31} & \cdots \\ \vdots & \ddots & O & O & \vdots & \ddots & \vdots & \vdots \\ \vdots & & \vdots & O & {t}_{\mathrm{31,3}m} & \cdots & {t}_{3m,3m}+{v}_{3m,3m}^{l} & \vdots \\ O & \cdots & O & \vdots & \vdots & \cdots & \cdots & \ddots \end{array}\right).\end{eqnarray} \tag{ 22 }$$

Since it is represented based on FE-DVR, the structure of Eq. (8) is also like the potential-energy matrix,

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{\chi }_{i,m}}{r}\right|{g}_{l-1}{rE}(t)\left|\displaystyle \frac{{\chi }_{i^{\prime},m^{\prime} }}{r}\right\rangle ={b}_{{im},{im}}^{l},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\begin{array}{rcl}{h}_{l,l-1}(t) & = & \left\langle \displaystyle \frac{{X}^{{\rm{T}}}(r)}{r}\right|{g}_{l-1}{rE}(t)\left|\displaystyle \frac{X(r)}{r}\right\rangle \\ & = & \left(\begin{array}{cccc}{b}_{11,11}^{l} & 0 & \cdots & 0\\ 0 & {b}_{12,12}^{l} & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & {b}_{{im},{im}}^{l}\end{array}\right).\end{array}\end{eqnarray} \tag{ 24 }$$

For ${h}_{l,l-1}(t)={h}_{l-1,l}(t)$, ${b}_{l}(t)$ in B(t) is expanded into a sparse symmetric matrix and expressed as

$$\begin{eqnarray}{b}_{l}(t)=\left(\begin{array}{ccccccc}0 & \cdots & 0 & {b}_{11,11}^{l} & 0 & \cdots & 0\\ \vdots & \ddots & \vdots & 0 & {b}_{12,12}^{l} & & \vdots \\ 0 & \cdots & 0 & \vdots & & \ddots & 0\\ {b}_{11,11}^{l} & 0 & \cdots & 0 & 0 & \cdots & {b}_{\mathrm{im},\mathrm{im}}^{l}\\ 0 & {b}_{12,12}^{l} & & \vdots & \ddots & & 0\\ \vdots & & \ddots & 0 & & \ddots & \vdots \\ 0 & \cdots & 0 & {b}_{\mathrm{im},\mathrm{im}}^{l} & 0 & \cdots & 0\end{array}\right).\end{eqnarray} \tag{ 25 }$$

Obviously, ${b}_{l}(t)$ is an ideal matrix for the SIL propagator in the perspective of efficiency. The number of nonzero elements in ${b}_{l}(t)$ is equal to the size of diagonal element of ${b}_{l}(t)$. This suggests that the ideal method of propagating partial waves in Eq. (17) is the SIL propagator. For the SIL propagator, an orthonormal Krylov subspace is constructed through recursive relations for each time step,

$$\begin{eqnarray}&&\left|{p}_{0}\right\rangle ={\Phi }_{l}(t)=\left(\begin{array}{l}{C}_{l-1}^{{\rm{T}}}(t)\\ {C}_{l}^{{\rm{T}}}(t)\end{array}\right),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\beta }_{0}\left|{p}_{1}\right\rangle =H\left|{p}_{0}\right\rangle -{\alpha }_{0}\left|{p}_{0}\right\rangle,\\ {\beta }_{j}\left|{p}_{j+1}\right\rangle =H\left|{p}_{j}\right\rangle -{\alpha }_{j}\left|{p}_{j}\right\rangle -{\beta }_{j-1}\left|{p}_{j-1}\right\rangle,\end{array}\end{eqnarray} \tag{ 27 }$$

where $H={b}_{l}(t)$, ${\alpha }_{j}$ and ${\beta }_{j}$ are expressed as

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{j}=\left\langle {p}_{j}\right|H\left|{p}_{j}\right\rangle,\\ {\beta }_{j}=\left\langle {p}_{j+1}\right|H\left|{p}_{j}\right\rangle .\end{array}\end{eqnarray} \tag{ 28 }$$

In the subspace ${P}_{N}=\left\{\left|{p}_{0}\right\rangle,\left|{p}_{1}\right\rangle,\ldots,\left|{p}_{N-1}\right\rangle \right\}$, ${H}_{{ij}}^{P}=\left\langle {p}_{i}\right|H\left|{p}_{j}\right\rangle$. N is the size of P_N. According to relations in Eqs. (27) and (28), H^p is a real tridiagonal matrix. The partial waves are propagating as follows:

$$\begin{eqnarray}&&{\Phi }_{l}(t+\delta t)\approx {U}^{P}{\Phi }_{l}(t),\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&{U}^{p}=\exp (-{\rm{i}}\delta {{tH}}^{p})=\mathop{\sum }\limits_{N+1}\exp (-{\rm{i}}\delta {{th}}_{i})\left|{Z}_{i}\right\rangle \left\langle {Z}_{i}\right|,\end{eqnarray} \tag{ 30 }$$

where $\left|{Z}_{i}\right\rangle$ and h_i are the corresponding eigenvalue and eigenvector of ${H}^{p}$. Because it is mainly affected by the operator of ${b}_{l}(t)$ times the vector p_j in Krylov subspace P_N, the calculation scale of the SIL operator is proportional to

$$\begin{eqnarray}&&{S}_{B}=N\times ({n}_{{\rm{e}}}{n}_{{\rm{g}}}-{n}_{{\rm{g}}}-1).\end{eqnarray} \tag{ 31 }$$

It is noted that for propagating the partial waves in Eq. (16) in the SIL frame, the calculation scale is proportional to

$$\begin{eqnarray}&&{S}_{A}=N\times ({n}_{{\rm{e}}}({n}_{{\rm{g}}}{n}_{{\rm{g}}}-1)-1),\end{eqnarray} \tag{ 32 }$$

because the number of nonzero elements in ${h}_{l,l}$ is ${n}_{{\rm{e}}}({n}_{{\rm{g}}}{n}_{{\rm{g}}}-1)-1$. Comparing with the ideal matrix like ${b}_{l}(t)$, the calculation scale of the SIL propagator for ${h}_{l,l}$ increases n_g times approximately when n_e is large. It indicates that the SIL propagator is not a good choice in this case. To overcome this difficulty, in the numerical calculation we decompose the matrix ${h}_{l,l}$ into smaller sub-blocks, and then the sub-blocks are diagonalized as given in Ref. [26]. Compared with diagonalizing the whole matrix of ${h}_{l,l}$ directly, this method can greatly reduce the amount of computation. Meanwhile, the matrix ${h}_{l,l}$ is time-independent, and the operations of decomposition and diagonalization of ${h}_{l,l}$ only need performing once before time-dependent propagating, which further reduces the amount of computation.

In the following, we mainly apply the present solver to the 3D TDSE of HHG from atomic hydrogen in a strong laser field. The parallel program is implemented in the FORTRAN language with openMP. The calculations are carried out on a workstation with two processors (Intel Xeon CPU, E5-2696v4, 22 cores, 2.20 GHz). In the calculations, the laser field is 800 nm in wavelength and 1.0×10¹4 W/cm² in power and has an envelope of ${\cos }^{2}(\pi t/{nT})$, where T is the time of an optical cycle, and n = 10, the full width of half maximum of the laser pulse is about 13.34 fs. The initial state of hydrogen atom is 1 s. The parameter of time is set to be T_total = 1200.0 a.u. and ${\rm{\Delta }}t=0.01$ a.u. The radial space region is in [0,200.0] a.u.

As is well known, the reliability of the solver is very important for simulating the real physical process of atoms interacted with strong laser fields. Figure 1 shows the HHG obtained from the present program, compared with the result of well-established SCID-TDSE,^[18] which is a program of time-dependent solution of 1-electron atomic Schrödinger equation in a strong laser field, and performed as well as LZH-DICP.^[24] The high-harmonic spectrum of SCID-TDSE is calculated from the Fourier transform of the real part of the expectation value of the dipole response. The high harmonic spectrum of our result is calculated from the Fourier transform of the real part of the expectation value of the time-dependent induced dipole acceleration response.^[34] In this calculation, the parameters are chosen to be n_g = 6, n_e = 200, L_max = 100, and N = 3. For obtaining a convergent result with the program SCID-TDSE, the parameters are ${\rm{\Delta }}r=0.3$ a.u., L_max = 100, and ${\rm{\Delta }}t=0.0005$ a.u. The consistency between the two results indicates that our solver is reliable. The visible difference between the two high-harmonic spectra in the high energy region is due to the fact that the precision of HHG with high energy obtained from the dipole response is not enough. It is noted that our program is several times faster than the SCID-TDSE when the two programs run on the same workstation with 32 threads, which is in line with expectations for the reason why the FDM of the program SCID-TDSE requires smaller time steps to ensure the convergence.^[35] The present solver can also be used to obtain the HHG of Ar atom in the same laser field with a model potential.^[36] There is a need for more FEs (n_e = 320) in the radial space for the steeper potential of Ar atom near nucleus. The power of the HHG of Ar atom as shown in Fig. 1(b) is greater than that of hydrogen atom for the larger tunnel rate of 3p state with higher angular momentum.^[37]

Zoom In Zoom Out Reset image size

The precision of the present solver is mainly affected by the discretization of space and time. The precision of the discretization of space based on the FE-DVR has been discussed in detail in Ref. [27]. The inappropriate values of n_e and n_g will increase the error of calculation as done by underfitting or overfitting. The error caused by the discretization of time has two main sources, one is from the operator splitting of exponential function in Eq. (15) which is proportional to ${\rm{\Delta }}{t}^{3}$, and the other is from the propagation based on the SIL propagator in Eqs. (29) and (30) that relies on the order of the Krylov subspace, and expressed as

$$\begin{eqnarray}&&\varepsilon \approx \displaystyle \frac{{\rm{\Delta }}{t}^{N-1}}{(N-1)!}\displaystyle \underset{i=0}{\overset{N-1}{\prod }}{\beta }_{i}.\end{eqnarray} \tag{ 33 }$$

To investigate the total error of the solver, the HHG of the calculations with different orders of the Krylov subspace is shown in Fig. 2(a). Figure 2(b) shows the comparison among the phases of HHG which are more sensitive to the accurate. The calculations are carried out with 32 threads. The parameters are n_g = 6, n_e = 200, and L_max = 200. Figure 2(a) and 2(b) show clearly that the results are convergent when $N\ge 3$, since the precision of the SIL propagator is consistent with the precision of the operator splitting of exponential function when N = 3 according to Eq. (31).

Zoom In Zoom Out Reset image size

It is noticeable that splitting out the centrifugal potential from the SIL propagator in our scheme is also one of the reasons of fast convergence with small order of Krylov subspace, which has the same result as that of Ref. [35]. The requirement for the small order of Krylov subspace means that a lot of time consumption of calculation can be saved according to Eq. (31). It is known that the requirement of accuracy increases with the increase of the laser intensity. So we give out the phase of HHG of Hydrogen atom in the laser field with its intensity up to 1.0×10¹⁵ W/cm² as shown in Fig. 3(a). It is obvious that three orders of Krylov subspace can meet the accuracy requirement. Three orders of Krylov subspace can also meet the accuracy requirement for Ar atom in an intense laser field as shown in Fig. 3(b).

Zoom In Zoom Out Reset image size

For a solver, propagating the wave packet efficiently with a dense matrix such as ${h}_{l,l}$ is a tricky problem. The propagator for the relatively dense matrix ${h}_{l,l}$ in our solver is not the SIL propagator, its time consumption does not dependent on the order of the Krylov subspace. It can be seen in Fig. 4(a) that the total time consumption of calculations varies with the order of the Krylov subspace as a linear function of N, specifically $f(N)=a\times N+b$. Thus, b = 42.27 (seconds) is the time consumption relating to the matrix ${h}_{l,l}$, and a × N (a = 11.63 (seconds)) is the time consumption for the matrix ${b}_{l}(t)$. If the SIL propagator is adopted for the matrix ${h}_{l,l}$, then based on Eqs. (31) and (32) it can be deduced that the time consumption can be given as $f(N)\approx (a/2)\times {n}_{{\rm{g}}}\times N$. For the cases of N = 3 and ${n}_{{\rm{g}}}=6$, the time consumption for the matrix ${h}_{l,l}$ is about 104.67 seconds, which is more time consuming than that of our method for the matrix ${h}_{l,l}$. This can make a significant statement that the propagation for the matrix ${h}_{l,l}$ in our solver is high efficient.

Zoom In Zoom Out Reset image size

Finally, the time consumption varying with the number of threads in the calculations of HHG is studied as shown in Fig. 4(b) with the parameters being L_max = 200, n_e = 200, and n_g = 6. This figure shows almost linearly scaling up to 32 threads, which indicates that the parallelism of the present solver is valid.

We have described an efficient parallel solver of TDSE simulating the interaction of atoms with a linearly polarized laser field. In this solver, the electronic wave function is expanded in terms of FE-DVR and the spherical harmonics. According to the RSPF, we implement the program in parallel, and improve the efficiency by separating out the ideal sparse sub-matrices from the Hamiltonian matrix for the SIL propagator. By decomposing and diagonalizing the sub-matrices including kinetic operator representation before time-dependent propagating, we can further improve the computational efficiency. Comparing with the split-Lanczos propagator,^[35] the efficiency of our solver is high because the SIL propagator is applied to the sparser matrix ${b}_{l}(t)$ instead of the matrices including kinetic operator. We use the solver to simulate HHG of atomic hydrogen in an intense laser field. The results show that the solver of TDSE has less time consumption and better parallel efficiency. In addition, it is expected that the present solver combined with the Wigner rotation technique should be utilized for simulating the interaction of atoms with circularly polarized intense laser field. For the TDSE for a two-electron system in spherical coordinates, a lot of block sub-matrices exist in the Hamiltonian matrix as a result of transition limits or parity conservation,^[20,38,39] our adopted method of separating and constructing block diagonal matrices may also increase the efficiency of solving the TDSE of a two-electron system interaction with an intense laser field.

We thank the High Performance Computing Center of Jilin University for the supercomputer time.