Efficient Time-Dependent Method for Strong-Field Ionization of Atoms with Smoothly Varying Radial Steps

Abstract

We present an efficient numerical method to solve the time-dependent Schrödinger equation in the single-active electron picture for atoms interacting with intense optical laser fields. Our approach is based on a non-uniform radial grid with smoothly increasing steps for the electron distance from the residual ion. We study the accuracy and efficiency of the method, as well as its applicability to investigate strong-field ionization phenomena, the process of high-order harmonic generation, and the dynamics of highly excited Rydberg states.

1. Introduction

The development of intense and tunable optical light sources has revolutionized atomic physics by revealing new strong-field phenomena, e.g., high-order harmonic generation (HHG) [1], above-threshold ionization (ATI) [2], Rydberg- state stabilization [3,4], or the transfer of the light’s orbital angular momentum (OAM) via optical vortex beams [5,6]. Even in atomic systems such as hydrogen, helium, the alkalis, and to some extent the heavy noble gases, analytical models to describe strong-field processes require the introduction of important approximations. Although numerous numerical methods, with various degrees of sophistication, have been designed in recent decades to treat atoms in optical fields, computing the electron dynamics in high-intensity long-wavelength pulses remains a challenge, even in the single-active electron (SAE) approach.

Recently, time-dependent spectral techniques [7,8,9], based on functions with compact support, e.g., B-splines [8,10,11,12] or Finite Element with Discrete Variable Representation (FE-DVR) [13,14], have been developed to treat multi-electron systems. Despite their impressive capabilities, these approaches remain computationally demanding and non-trivial to implement. In addition, they can hardly describe high-intensity optical pulses, besides a few exceptions, such as the R-matrix with time dependence (RMT) [7], or methods based on infinite exterior complex scaling (iECS) [15]. These methods also require handling large-size matrices, which can lead to additional complications and limitations.

On the other hand, SAE approaches based on finite-difference methods remain among the most favorite methods to implement, while enabling an accurate treatment of a plethora of strong-field processes where the single-electron dynamics dominates. Indeed, the Crank–Nicholson method [16], combined with a three-point formula to evaluate the action of the kinetic energy operator, allows one to produce highly efficient propagation schemes. Moreover, these methods can be combined with sophisticated techniques, such as flux methods [17,18], to produce some of the most efficient time-dependent methods.

An important limiting aspect of methods based on finite differences lies in the choice of the radial grid. Typically, a uniform grid with a fine radial mesh is needed to describe the highly localized electronic wave function at short distances, where the electron possesses substantial kinetic energy. Otherwise, bound-state energies can rapidly become inaccurate. At large radial distances, however, low-energy electrons, as produced in many strong-field ionization processes with optical pulses [19,20,21], do not require a high density of radial points due to their relatively long wavelength. The same principle applies in HHG, since the recolliding electron responsible for the light emission in the high-harmonic plateau region acquires its highest kinetic energy near the ionic core. Therefore, the need for a small radial mesh at short distances often imposes an unnecessary constraint at large distances if a uniform grid is employed. This restriction, however, does not exist with the use of compact functions since their density as a function of radial distance can easily be varied.

In this work, we present the implementation of a non-uniform radial grid with a finite-difference method to solve the time-dependent Schrödinger equation (TDSE) for atoms in optical fields. We build upon past methods developed for one-dimensional time-independent problems [22,23] to the three-dimensional time-dependent case, focusing on the study of atoms interacting with high-intensity optical pulses. We demonstrate considerable improvements in time propagation while at the same time conserving high accuracy. The validity of the scheme is demonstrated through benchmark calculations performed on the hydrogen atom interacting with infrared light. However, the method can be employed for any atom where the single-electron dynamics can be approximately described by a mean-field potential, such as the alkalis, helium, and often even the heavier noble gases. In particular, we show applications of the non-uniform grid to efficiently treat strong-field ionization, the computation of HHG spectra, and the dynamics of highly excited Rydberg states.

This manuscript is constructed as follows. In the next section, we describe the approach to solve the TDSE with a non-uniform radial grid, and we also discuss the theoretical framework needed to compute physical observables. In Section 3, we demonstrate the efficiency of the method and its potential applications in strong-field physics. We finish with a summary and our conclusions in Section 4.

Unless stated otherwise, atomic units are used throughout this manuscript.

2. Theoretical Approach

2.1. General Formalism Using a Uniform Radial Grid

In the SAE approach, the TDSE for an atomic system interacting with an electromagnetic field takes the form

$$-i{\partial }_t|\mathsf{\Psi }\left(t\right)\rangle =\widehat{H}\left(t\right)|\mathsf{\Psi }\left(t\right)\rangle ,$$

(1)

where $\widehat{H}\left(t\right)={\widehat{H}}_0+{\widehat{H}}_{\mathrm{int}}\left(t\right)$ is the time-dependent Hamiltonian of the system comprising the field-free Hamiltonian ${\widehat{H}}_0=\widehat{T}+\widehat{V}$, where $\widehat{T}$ and $\widehat{V}$ are the kinetic energy and the spherical electronic potential operators, respectively, while ${\widehat{H}}_{\mathrm{int}}\left(t\right)$ is the atom–field interaction Hamiltonian. Employing the dipole approximation and the Coulomb gauge, we express the interaction operator either in the length gauge (LG), ${\widehat{H}}_{\mathrm{int}}^L\left(t\right)=\mathit{r}·\mathit{E}\left(t\right)$, or in the velocity gauge (VG), ${\widehat{H}}_{\mathrm{int}}^V\left(t\right)=\mathit{p}·\mathit{A}\left(t\right)$, where $\mathit{A}\left(t\right)$ and $\mathit{E}\left(t\right)=-d\mathit{A}\left(t\right)/dt$ are the vector potential and electric field, respectively, while $\mathit{r}$ and $\mathit{p}$ are the electron position and momentum vectors. To simplify the presentation of the method in the following development, we consider linearly polarized light with the electric-field vector along the z-axis and define the vector potential as

$$\mathit{A}\left(t\right)=A_0{sin}^2\left(\frac{\omega t}{2N}\right)sin\left(\omega t\right)\widehat{z}$$

(2)

for $0\le t\le t_{\mathrm{end}}$, where $t_{\mathrm{end}}=2\pi N/\omega$ and N is the number of cycles. The method can be extended to encompass light with arbitrary polarization in a straightforward way.

The temporal evolution of the wave function from time $t_n=n\Delta t$ to $t_{n+1}=t_n+\Delta t$ is given by the relation $|\mathsf{\Psi }\left(t_{n+1}\right)\rangle =\widehat{U}(t_{n+1},t_n)|\mathsf{\Psi }\left(t_n\right)\rangle$, where $\widehat{U}(t_{n+1},t_n)$ is the time-evolution operator in the Schrödinger picture between $t_n$ and $t_{n+1}$. The exponential time propagator is further split into field-free and atom–field interaction components as

$$\widehat{U}(t_{n+1},t_n)\approx e^{-i{\widehat{H}}_0\frac{\Delta t}{2}}{e}^{-i{\widehat{H}}_{\mathrm{int}}(t_n+\Delta t/2)\Delta t}{e}^{-i{\widehat{H}}_0\frac{\Delta t}{2}},$$

(3)

resulting in a third-order method in the time step.

The action of these operators is evaluated by expanding the wave function in a truncated series of partial waves as

$$\mathsf{\Psi }(\mathit{r},t)=\frac{1}{r}\sum _{\ell =\left|m\right|}^{{\ell }_{\max }}{u}_{\ell }(r,t)Y_{\ell }^m\left({\mathsf{\Omega }}_e\right),$$

(4)

where r is the electron’s radial distance, ${\mathsf{\Omega }}_e$ denotes its spherical coordinates, and the $Y_{\ell }^m\left({\mathsf{\Omega }}_e\right)=\langle {\mathsf{\Omega }}_e|\ell m\rangle$ are the spherical harmonics with ℓ and m representing the electron’s angular momentum and its projection on the z-axis, respectively. In Equation (4), we impose the boundary condition $u_{\ell }(0,t)=0$ at any time t, and m is fixed by the magnetic quantum number of the initial state.

The field-free Hamiltonian in spherical coordinates is represented by

$${\widehat{H}}_0=-\frac{1}{2}\left[\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial }{\partial r}\right)-\frac{{\widehat{\ell }}^2}{r^2}\right]+V\left(r\right),$$

(5)

where $\widehat{\ell }$ is the one-particle orbital angular momentum operator acting on the angular part of the wave function. The action of ${\widehat{H}}_0$ on $\left|\mathsf{\Psi }\right(t)\rangle$, expressed in the $|r\ell m\rangle$ basis, is given by

$$\langle r\ell m|{\widehat{H}}_0|\mathsf{\Psi }\left(t\right)\rangle =-\frac{1}{2}\frac{{\partial }^2{u}_{\ell }(r,t)}{\partial r^2}+V_{\mathrm{eff}}\left(r\right)u_{\ell }(r,t),$$

(6)

where $V_{\mathrm{eff}}\left(r\right)=V\left(r\right)+\ell (\ell +1)/2r^2$.

We initially introduce a standard uniform radial grid $r_j=(j-1)h$ of constant step h for $j=1,\cdots ,M$, and the vectors ${\mathit{u}}_{\ell }\left(t\right)$ with components $u_{\ell ,j}\left(t\right)=u_{\ell }(r_j,t)$. Using the second-order central difference formula to approximate the radial kinetic energy term for equidistant grid points gives

$$\frac{{\partial }^2{u}_{\ell }(r_j,t)}{\partial r^2}=\frac{u_{\ell ,j+1}\left(t\right)-2u_{\ell ,j}\left(t\right)+u_{\ell ,j-1}\left(t\right)}{h^2}+\mathcal{O}\left(h^2\right).$$

(7)

This leads to the following non-zero elements of the field-free finite-difference Hamiltonian matrix

$${\left(H_0\right)}_{j,j^{\prime }}^{\ell ,{\ell }^{\prime }}={\delta }_{\ell ,{\ell }^{\prime }}\left\{\begin{array}{c}\hfill {\displaystyle -\frac{1}{2h^2},j^{\prime }=j\pm 1,}\hfill \\ \hfill {\displaystyle \frac{1}{h^2}+V_{\mathrm{eff}}\left(r_j\right),j^{\prime }=j,}\hfill \end{array}\right.$$

(8)

which exhibits a block-diagonal structure in the partial waves $\ell ,{\ell }^{\prime }$ and a tridiagonal form in the radial grid points. Utilizing Cayley’s transform [24] of the complex exponential time-evolution operator for the field-free Hamiltonian, i.e.,

$$e^{-i{\widehat{H}}_0\frac{\Delta t}{2}}\approx {[\widehat{I}+i{\widehat{H}}_0\Delta t/4]}^{-1}[\widehat{I}-i{\widehat{H}}_0\Delta t/4]+\mathcal{O}(\Delta t^3),$$

(9)

provides an unconditionally stable third-order method whose action can be evaluated by either inverting a tridiagonal matrix or by solving a tridiagonal system of linear equations using the Thomas algorithm.

Moving to the atom–field interaction Hamiltonian in the length gauge, ${\widehat{H}}_{\mathrm{int}}^L\left(t\right)=E\left(t\right)rcos\theta$, its action on $\left|\mathsf{\Psi }\right(t)\rangle$ is expressed as

$$\langle r\ell m|{\widehat{H}}_{\mathrm{int}}^L\left(t\right)|\mathsf{\Psi }\left(t\right)\rangle =E\left(t\right)r[c_{\ell }^m{u}_{\ell +1}(r,t)+c_{\ell -1}^m{u}_{\ell -1}(r,t)],$$

(10)

where the coefficients $c_{\ell }^m$ are given by

$$c_{\ell }^m={\left[\frac{(\ell -m+1)(\ell +m+1)}{(2\ell +1)(2\ell +3)}\right]}^{\frac{1}{2}}.$$

(11)

The non-zero elements of the interaction Hamiltonian take the form

$${\left(H_{\mathrm{int}}^L\right)}_{j,j^{\prime }}^{\ell ,{\ell }^{\prime }}={\delta }_{j,j^{\prime }}E\left(t\right)r_j\left\{\begin{array}{c}\hfill {\displaystyle c_{\ell }^m,{\ell }^{\prime }=\ell +1,}\hfill \\ \hfill c_{\ell -1}^m,{\ell }^{\prime }=\ell -1,\hfill \end{array}\right.$$

(12)

exhibiting a block-diagonal structure in the radial coordinate and a tridiagonal form in the angular momentum quantum number. The same approach used in the case of the field-free Hamiltonian can thus be employed to evaluate its action on the wave function.

The form of the velocity–gauge interaction Hamiltonian in spherical coordinates is given by

$${\widehat{H}}_{\mathrm{int}}^V\left(t\right)=-iA\left(t\right)\left(cos\theta \frac{\partial }{\partial r}-\frac{sin\theta }{r}\frac{\partial }{\partial \theta }\right).$$

(13)

Using a central difference formula to evaluate the first derivative, the non-zero elements of the finite difference interaction Hamiltonian become

$${\left(H_{\mathrm{int}}^V\right)}_{j,j^{\prime }}^{\ell ,{\ell }^{\prime }}=-iA\left(t\right)\left\{\begin{array}{c}{\displaystyle \pm \frac{c_{\ell }^m}{2h},{\ell }^{\prime }=\ell +1,j=j^{\prime }\pm 1,}\hfill \\ {\displaystyle (\ell +2)\frac{c_{\ell }^m}{r_j},{\ell }^{\prime }=\ell +1,j=j^{\prime },}\hfill \\ {\displaystyle \pm \frac{c_{\ell -1}^m}{2h},{\ell }^{\prime }=\ell -1,j=j^{\prime }\pm 1,}\hfill \\ {\displaystyle -(\ell -1)\frac{c_{\ell -1}^m}{r_j},{\ell }^{\prime }=\ell -1,j=j^{\prime }.}\hfill \end{array}\right.$$

(14)

Due to the radial derivative, the velocity gauge requires more computational operations than the length gauge. However, since fewer partial waves are often sufficient in the velocity gauge to achieve convergence, this gauge becomes more efficient at high intensity and long wavelengths [25,26].

Finally, a maximum grid size $R_{\max }$ should be chosen such that the propagated wave packet remains within the grid at all times. Since the description of the electron dynamics is only required in the region where the wave packet has expanded, it proves efficient to employ a dynamically increasing grid size $R_{\max }\left(t\right)$ that contains the wave packet at any time, but avoids unnecessary computation at distances where the wave function’s value is negligible. Here, we expand our grid as $R_{\max }(t+\Delta t)=R_{\max }\left(t\right)+\Delta R$, whenever it is found that the wave packet probability amplitude has increased above a certain accuracy parameter, $\left|\mathsf{\Psi }\right(\mathit{r},t\left)\right|\ge \epsilon$, in a region just below $R_{\max }\left(t\right)$, i.e., for $(R_{\max }\left(t\right)-\Delta R)<r<R_{\max }\left(t\right)$, where $\Delta R$ is a radial interval. As a result, one can still obtain highly converged results while saving a considerable amount of computational resources.

2.2. Modification of the Formalism with a Non-Uniform Radial Grid

2.2.1. Finite Difference with Smoothly Varying Radial Steps

The efficiency of the method outlined in Section 2.1 stems from the tridiagonal structure of the kinetic energy operator, originating from the three-point kinetic energy formula in (7). Therefore, it is desirable to seek a method preserving the three-point finite scheme while enabling flexibility in the choice of the radial steps.

To accomplish this goal, we introduce a non-uniform radial grid with non-equidistant steps $h_j=r_{j+1}-r_j$, for $j=1,\cdots ,M$, with $r_1=0$, as illustrated in Figure 1. Omitting the time dependence to simplify the expressions, we first expand $u_{\ell ,j+1}=u_{\ell }\left(r_{j+1}\right)$ and $u_{\ell ,j-1}=u_{\ell }\left(r_{j-1}\right)$ in Taylor series around $r_j$:

$$\begin{array}{ccc}\hfill u_{\ell ,j+1}& =& \hfill u_{\ell ,j}+h_j{u}_{\ell ,j}^{\prime }+\frac{h_j^2}{2}{u}_{\ell ,j}^{\prime \prime }+\frac{h_j^3}{6}{u}_{\ell ,j}^{\prime \prime \prime }+\frac{h_j^4}{24}{u}_{\ell ,j}^{\left(4\right)}+\mathcal{O}\left(h_j^5\right),\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill u_{\ell ,j-1}& =& \hfill u_{\ell ,j}-h_{j-1}{u}_{\ell ,j}^{\prime }+\frac{h_{j-1}^2}{2}{u}_{\ell ,j}^{\prime \prime }-\frac{h_{j-1}^3}{6}{u}_{\ell ,j}^{\prime \prime \prime }+\frac{h_{j-1}^4}{24}{u}_{\ell ,j}^{\left(4\right)}+\mathcal{O}\left(h_{j-1}^5\right),\hfill \end{array}$$

(16)

where the primed $u_{\ell ,j}$ represent derivatives of $u_{\ell }$ evaluated at $r_j$. We solve for $u_{\ell ,j}^{\prime \prime }$ by taking $h_{j-1}\times \left(15\right)+h_j\times \left(16\right)$. This leads to

$$\begin{array}{ccc}\hfill u_{\ell ,j}^{\prime \prime }& =& \hfill \frac{2\left(h_{j-1}{u}_{\ell ,j+1}-(h_j+h_{j-1})u_{\ell ,j}+h_j{u}_{\ell ,j-1}\right)}{h_j{h}_{j-1}^2+h_j^2{h}_{j-1}}\hfill \\ & & \hfill -\frac{(h_j-h_{j-1})}{3}{u}_{\ell ,j}^{\prime \prime \prime }+\mathcal{O}\left(h_j^2\right)+\mathcal{O}\left(h_{j-1}^2\right).\hfill \end{array}$$

(17)

Equation (17) gives the correct expression for $u_{\ell ,j}^{\prime \prime }$ to $\mathcal{O}(h_j-h_{j-1})$, i.e., to first order in the difference in the grid increment. For a uniform grid, $h_j=h_{j-1}$, and we recover the second-order central difference formula in (7), i.e.,

$$u_{\ell ,j}^{\prime \prime }=\frac{u_{\ell ,j+1}-2u_{\ell ,j}+u_{\ell ,j-1}}{h^2}+\mathcal{O}\left(h^2\right).$$

(18)

The subtle cancellation of first-order terms due to the symmetry of the uniform grid is thus responsible for the second-order accuracy of the method. Any increment of the form $h_j-h_{j-1}=\mathcal{O}\left(h_{j-1}\right)$, even if applied at just one grid point to bridge grids with distinct step sizes, will introduce a significant error that will propagate across space and time, ultimately undermining the accuracy of the method.

Consequently, it is essential to choose a grid increment such that $h_j-h_{j-1}=\mathcal{O}\left(h_{j-1}^2\right)$ to maintain second-order accuracy. Following the idea introduced in [22], we employ grid steps defined by the recursive relation

$$h_j=h_{j-1}(1+\alpha h_{j-1}),$$

(19)

where $\alpha$ is a fixed parameter controlling the rate of incremental change in the radial step size. At larger distances, we impose an upper limit $h_{\max }$ on the radial step size to be able to describe electrons with a maximum momentum $p_0$. Consequently, beyond a grid index $j_0$, we use again a standard uniform grid where $h_j=h_{\max }$ for all $j\ge j_0$ (see Figure 1). Typically, about ten radial points should span the electronic wavelength for an accurate representation of the wave function, leading to $h_{\max }\approx \pi /5p_0$. For electron energies up to about 20 eV, this corresponds to a maximum grid size of $h_{\max }\approx 0.5$ a.u. This value should be contrasted with the radial step of $h\approx 0.01$ a.u. needed to reproduce the hydrogen ground-state energy with high precision.

The non-zero elements of the field-free Hamiltonian for the general case of a non-uniform grid become

$${\left(H_0\right)}_{j,j^{\prime }}^{\ell ,{\ell }^{\prime }}={\delta }_{\ell ,{\ell }^{\prime }}\left\{\begin{array}{c}{\displaystyle -\frac{1}{2}\frac{1}{L_j^2}\frac{1}{h_j},j^{\prime }=j+1,}\hfill \\ {\displaystyle -\frac{1}{2}\frac{1}{L_j^2}\frac{1}{h_{j-1}},j^{\prime }=j-1,}\hfill \\ {\displaystyle \frac{1}{2}\frac{1}{L_j^2}\left(\frac{1}{h_{j-1}}+\frac{1}{h_j}\right)+V_{\mathrm{eff}}\left(r_j\right),j^{\prime }=j,}\hfill \end{array}\right.$$

(20)

where we introduced $L_j^2=(h_{j-1}+h_j)/2$. Although ${\widehat{H}}_0$ keeps its tridiagonal structure, it is no longer symmetric. As described in the next section, this tridiagonal matrix is similar to a hermitian matrix through a diagonal change of the basis matrix. Hence, its eigenvalues are real, and the general properties resulting from the hermicity of the operators can be recovered by using the appropriate metric space (see Section 2.2.2).

The interaction operator in the length gauge (12) is unaffected by the change to a non-uniform grid. However, the velocity–gauge operator requires the evaluation of the first radial derivative. Repeating our previous approach, we can isolate $u_{\ell ,j}^{\prime }$ by taking $h_{j-1}^2\times \left(15\right)-h_j^2\times \left(16\right)$, leading to a three-point formula [22] with an error of order $\mathcal{O}\left(h_j{h}_{j-1}\right)$. In our case, the application of a three-point formula to evaluate the first derivative is unnecessary since the same accuracy can already be achieved with a simpler two-point central-difference formula. Indeed, subtracting directly (15) and (16), we obtain

$$\begin{array}{ccc}\hfill u_{\ell ,j}^{\prime }& =& \frac{u_{\ell ,j+1}-u_{\ell ,j-1}}{h_{j-1}+h_j}+\frac{h_j^2-h_{j-1}^2}{2}{u}_{\ell ,j}^{\prime \prime }+\mathcal{O}\left(h_{j-1}^2\right)\hfill \\ & & +\mathcal{O}\left(h_j^2\right)+\mathcal{O}\left(h_{j-1}{h}_j\right)\hfill \end{array}$$

(21)

for the general case. With our choice of grid in (19), $h_j^2-h_{j-1}^2=\mathcal{O}\left(h_{j-1}^3\right)$, and Equation (21) is thus also accurate to second order. Using the latter formula, the non-zero elements of the interaction Hamiltonian in the velocity gauge become

$${\left(H_{\mathrm{int}}^V\right)}_{j,j^{\prime }}^{\ell ,{\ell }^{\prime }}=-iA\left(t\right)\left\{\begin{array}{c}{\displaystyle \pm \frac{1}{2}\frac{c_{\ell }^m}{L_j},{\ell }^{\prime }=\ell +1,j=j^{\prime }\pm 1,}\hfill \\ {\displaystyle (\ell +2)\frac{c_{\ell }^m}{r_j},{\ell }^{\prime }=\ell +1,j=j^{\prime },}\hfill \\ {\displaystyle \pm \frac{1}{2}\frac{c_{\ell -1}^m}{L_j},{\ell }^{\prime }=\ell -1,j=j^{\prime }\pm 1,}\hfill \\ {\displaystyle -(\ell -1)\frac{c_{\ell -1}^m}{r_j},{\ell }^{\prime }=\ell -1,j=j^{\prime }.}\hfill \end{array}\right.$$

(22)

The implementation of the non-uniform grid is seen to involve straightforward modifications of the standard SAE method, but it can yield a considerable improvement in computational efficiency. The method, however, necessitates the utilization of a consistent inner product rule to normalize eigenstates and to obtain physical observables, as detailed below.

2.2.2. Modified Inner Product Rule with the Non-Uniform Grid

Let us consider the finite-difference field-free Hamiltonian ${\widehat{H}}_0$ in the case of a non-uniform radial grid (20). Since the Hamiltonian is tridiagonal, and the upper and lower diagonals have the same sign, ${\widehat{H}}_0$ can be diagonalized, ${\widehat{H}}_0|{\mathsf{\Psi }}_{n,\ell }\rangle =E_{n,\ell }|{\mathsf{\Psi }}_{n,\ell }\rangle$, with real eigenvalues $E_{n,\ell }$ and eigenstates $|{\mathsf{\Psi }}_{n,\ell }\rangle$. Here, we enforce fixed-boundary conditions on the wave function, i.e., $u_{\ell ,1}=u_{\ell ,M}=0$, at $r_1=0$ and at the maximum grid size $r_M$. The eigenstates $|{\mathsf{\Psi }}_{n,\ell }\rangle$ with $E_{n,\ell }<0$ correspond to the bound states of the atomic system.

As previously noted, ${\widehat{H}}_0$ is not hermitian for a non-uniform grid. It can, however, be transformed into a hermitian matrix ${\overline{H}}_0$ by a similar transformation using a diagonal matrix $\widehat{L}$ with elements $L_j$ for $j=1,\cdots ,M$ in each ℓ block. Specifically, ${\overline{H}}_0=\widehat{L}{\widehat{H}}_0{\widehat{L}}^{-1}$ with

$${\left(\widehat{L}\right)}_{j,j^{\prime }}^{\ell ,{\ell }^{\prime }}={\delta }_{\ell ,{\ell }^{\prime }}{\delta }_{j,j^{\prime }}{L}_j.$$

(23)

Introducing the functions $|{\overline{\mathsf{\Psi }}}_{n,\ell }\rangle =\widehat{L}|{\mathsf{\Psi }}_{n,\ell }\rangle$, $|{\overline{\mathsf{\Psi }}}_{n,\ell }\rangle$ is an eigenstate of ${\overline{H}}_0$ with eigenvalue $E_{n,\ell }$, since

$${\overline{H}}_0|{\overline{\mathsf{\Psi }}}_{n,\ell }\rangle =\widehat{L}{\widehat{H}}_0{\widehat{L}}^{-1}\widehat{L}|{\mathsf{\Psi }}_{n,\ell }\rangle =E_{n,\ell }|{\overline{\mathsf{\Psi }}}_{n,\ell }\rangle .$$

(24)

Hence, the eigenstates $|{\mathsf{\Psi }}_{n,\ell }\rangle$ of ${\overline{H}}_0$ can be normalized to form an orthonormal basis in the sense that

$$\langle {\overline{\mathsf{\Psi }}}_{m,\ell }|{\overline{\mathsf{\Psi }}}_{n,{\ell }^{\prime }}\rangle =\langle {\mathsf{\Psi }}_{m,\ell }|\widehat{M}|{\mathsf{\Psi }}_{n,{\ell }^{\prime }}\rangle ={\delta }_{m,n}{\delta }_{\ell ,{\ell }^{\prime }},$$

(25)

where $\widehat{M}={\widehat{L}}^2$ is a diagonal weighting matrix. Considering two given wave functions $|\mathsf{\Lambda }\rangle$ and $|\mathsf{\Omega }\rangle$, respectively, defined by their discrete radial components $u_{j,\ell }$ and $v_{j,\ell }$, their inner product is given by

$$\langle \overline{\mathsf{\Lambda }}|\overline{\mathsf{\Omega }}\rangle =\langle \mathsf{\Lambda }|\widehat{M}|\mathsf{\Omega }\rangle =\sum _{\ell =\left|m\right|}^{{\ell }_{\max }}\sum _{j=2}^{M-1}{L}_j^2{u}_{j,\ell }^{*}{v}_{j,\ell }.$$

(26)

In particular, this metric preserves the norm of the wave function through time propagation. Indeed, since ${\overline{H}}_{\mathrm{int}}^L=H_{\mathrm{int}}^L$ and ${\overline{H}}_{\mathrm{int}}^V=\widehat{L}{H}_{\mathrm{int}}^V{\widehat{L}}^{-1}$ are both hermitian, the time-dependent Hamiltonian $\overline{H}\left(t\right)=\widehat{L}\widehat{H}\left(t\right){\widehat{L}}^{-1}$ is hermitian in both gauges. If $\left|\mathsf{\Psi }\right(t)\rangle$ is a solution of Equation (1), then $|\overline{\mathsf{\Psi }}\left(t\right)\rangle =\widehat{L}|\mathsf{\Psi }\left(t\right)\rangle$ is a solution of $-i{\partial }_t|\overline{\mathsf{\Psi }}\left(t\right)\rangle =\overline{H}\left(t\right)|\overline{\mathsf{\Psi }}\left(t\right)\rangle$, and the norm of $\left|\mathsf{\Psi }\right(t)\rangle$ is thus conserved under the inner product rule (26).

2.2.3. Calculation of Physical Observables

The modification of the inner product rule necessitates some adjustments when evaluating physical observables. Let us consider the time-dependent wave function $\left|\mathsf{\Psi }\right(t)\rangle$ solution of the TDSE in Equation (1) along with an observable of the form $\widehat{\mathcal{O}}(\nabla ,{\nabla }^2,{\mathit{r}}^n)$, whose action is evaluated on the non-uniform grid using the finite-difference formula introduced in Section 2.2. This observable can be transformed into a hermitian operator through the similarity transformation, $\overline{\mathcal{O}}=\widehat{L}\mathcal{O}{\widehat{L}}^{-1}$. Its expectation value is calculated as

$$\langle \overline{\mathsf{\Psi }}\left(t\right)|\overline{\mathcal{O}}|\overline{\mathsf{\Psi }}\left(t\right)\rangle =\langle \mathsf{\Psi }\left(t\right)|\widehat{M}\widehat{\mathcal{O}}|\mathsf{\Psi }\left(t\right)\rangle .$$

(27)

Moreover, the probability of finding the electron at the end of the pulse in an excited atomic bound state $|{\mathsf{\Psi }}_{n,\ell }\rangle$ is given by $|\langle {\mathsf{\Psi }}_{n,\ell }|\widehat{M}{|\mathsf{\Psi }\left(t\right)\rangle |}^2$. To calculate physical quantities related to asymptotically free states of the electron, one expands the eigenfunction of the atomic Hamiltonian, corresponding to a photoelectron with asymptotic momentum $\mathit{k}$, in partial waves as

$${\mathsf{\Phi }}_{\mathit{k}}^{(-)}=\frac{4\pi }{r}\sum _{\ell m}{i}^{\ell }{e}^{-i{\delta }_{E\ell }}{w}_{E\ell }\left(r\right)Y_{\ell }^{m\ast }\left({\mathsf{\Omega }}_{\mathit{k}}\right)Y_{\ell }^m\left({\mathsf{\Omega }}_e\right).$$

(28)

In the above equation, ${\mathsf{\Omega }}_{\mathit{k}}$ are the spherical angles that characterize the direction of $\mathit{k}$, ${\delta }_{E\ell }={\delta }_{\ell }^0+{\sigma }_{\ell }$ is the scattering phase, including the short-range potential-scattering phase shift ${\delta }_{\ell }^0$ and the Coulomb phase ${\sigma }_{\ell }=\mathsf{\Gamma }(\ell +1-i/k)$. The radial functions $w_{E\ell }\left(r\right)$ are solutions of the time-independent equation

$$\frac{d^2{w}_{E\ell }\left(r\right)}{d{r}^2}+2\left[E-V_{\mathrm{eff}}\left(r\right)\right]w_{E\ell }\left(r\right)=0.$$

(29)

Here, it is essential to use (17) to evaluate the second derivative in (29). Introducing $w_{j,E\ell }=w_{E\ell }\left(r_j\right)$, we obtain the relation

$$\begin{array}{c}\hfill w_{j+1,E\ell }=\left(1+\frac{h_j}{h_{j-1}}\right)w_{j,E\ell }-\frac{h_j}{h_{j-1}}{w}_{j-1,E\ell }-2h_{j-1}{L}_j(E-V_{\mathrm{eff}}\left(r_j\right))w_{j,E\ell }\end{array}$$

(30)

for $j\ge 2$. This can be solved recursively, after fixing $w_{1,E\ell }=0$ and $w_{2,E\ell }=w_0$, where $w_0$ acts as an overall scaling factor of the wave function. Using energy-normalized scattering functions, the energy-resolved angular distribution calculated after the end of the pulse, $t=t_f$, is given by

$$\begin{array}{ccc}\hfill \frac{d\mathcal{P}}{dEd{\mathsf{\Omega }}_{\mathit{k}}}& =& |\langle {\mathsf{\Phi }}_{\mathit{k}}^{(-)}|\widehat{M}|\mathsf{\Psi }\left(t_f\right){\rangle |}^2,\hfill \\ \hfill \frac{d\mathcal{P}}{dEd{\mathsf{\Omega }}_{\mathit{k}}}& =& \hfill \frac{{\mathcal{P}}_E}{4\pi }\left[1+\sum _{L>0}{\beta }_L\left(E\right)P_{\ell }(cos{\theta }_{\mathit{k}})\right],\hfill \end{array}$$

(31)

where $E={\mathit{k}}^2/2$, $P_{\ell }\left(x\right)$ is a standard Legendre polynomial, and ${\theta }_{\mathit{k}}$ is the polar angle characterizing the direction of $\mathit{k}$. The angle-integrated probability density can be expressed as

$${\mathcal{P}}_E=\int d{\mathsf{\Omega }}_{\mathit{k}}\frac{d\mathcal{P}}{dEd{\mathsf{\Omega }}_{\mathit{k}}}=\sum _{\ell =\left|m\right|}^{{\ell }_{\max }}{\left|Z_{E\ell }\right|}^2,$$

(32)

where the coefficients

$$Z_{E\ell }=\sum _{j=2}^{M-1}{L}_j^2{w}_{j,E\ell }{u}_{j,\ell }$$

(33)

are evaluated after the end of the pulse. The total ionization probability is given by $W=\int {\mathcal{P}}_{E}dE$, and the anisotropy parameters take the form

$${\beta }_L\left(E\right)=\frac{1}{{\mathcal{P}}_E}\sum _{\ell {\ell }^{\prime }}{i}^{{\ell }^{\prime }-\ell }{e}^{i({\delta }_{E\ell }-{\delta }_{E{\ell }^{\prime }})}{\nu }_{\ell {\ell }^{\prime }m}^L{Z}_{E{\ell }^{\prime }}^{*}{Z}_{E\ell },$$

(34)

where ${\nu }_{\ell {\ell }^{\prime }m}^L=\sqrt{(2\ell +1)(2{\ell }^{\prime }+1)}{C}_{\ell m,{\ell }^{\prime }-m}^{L0}$ with $C_{\ell m,{\ell }^{\prime }-m}^{L0}$ denoting a Clebsch–Gordan coefficient.

3. Results and Discussion

3.1. Benchmark Study on Strong-Field Ionization

We start by investigating the efficiency and accuracy of the non-uniform grid with our focus on strong-field ionization of atomic hydrogen. We employ a 400 nm pulse with $N=6$ cycles at a peak intensity $I={10}^{14}$ W/cm². Considering at first the LG, a highly converged spectrum was achieved for a uniform grid with $h=0.01$, ${\ell }_{\max }=25$, $dt=0.01$, and dynamical grid parameters $\epsilon ={10}^{-10}$ and $\Delta R=20$ (see Section 2.1). The calculated ionization spectrum, shown in Figure 2a, will serve as a benchmark for assessing the accuracy of subsequent results.

The calculation was then repeated for uniform grids with different radial steps $h=0.1,0.25,$ and $0.5$ while keeping all other parameters unchanged. For each uniform radial grid, the initial ground state of hydrogen was recomputed by diagonalizing the field-free Hamiltonian in (8). The obtained spectra are shown in Figure 2a and reveal significant discrepancies with our benchmark calculation at $h=0.01$. Noticeable differences are already apparent for a radial step as small as $h=0.1$. These discrepancies stem from the inaccurate description of the ground-state wave function, which can also lead to incorrect spectral features. These results underscore the necessity of employing a uniform grid with small steps to reproduce accurate spectra.

We proceeded to replicate the calculation employing different non-uniform radial grids, as defined in Equation (19), setting an increment parameter of $\alpha =0.05$ throughout the following development. In all calculations, we fixed the initial radial step $h_1=0.01$ to maintain a good representation of the ground state, and we used various maximum mesh sizes, $h_{\max }=0.01,0.1,0.25,$ and $0.5$, while again keeping all other parameters unchanged. Therefore, the grid with $h_{\max }=0.01$ corresponds to our benchmark calculation and will be referred to as such throughout our discussion. The ionization spectra with different non-uniform grids are depicted in Figure 2b and show excellent agreement with the benchmark calculation. Upon closer inspection, slight discrepancies become apparent in the high-energy part of the spectrum for $h_{\max }=0.5$. To assess the accuracy of the calculations, we evaluated the relative error as ${ϵ}_r\left(E\right)=|{\mathcal{P}}_E^{\mathrm{ng}}-{\mathcal{P}}_E^{\mathrm{be}}|/{\mathcal{P}}_E^{\mathrm{be}}$, where ${\mathcal{P}}_E^{\mathrm{ng}}$ and ${\mathcal{P}}_E^{\mathrm{be}}$ represent the energy-resolved ionization probability computed with the non-uniform grid and benchmark (uniform grid) calculations, respectively. In addition, we introduced the maximum relative error ${ϵ}_r^{\max }\left(\mathcal{I}\right)={max}_{\{E\in \mathcal{I}\}}{ϵ}_r\left(E\right)$ over an energy interval $\mathcal{I}$, and the relative error on the total ionization probability ${ϵ}_r^{\mathrm{ion}}=|W^{\mathrm{ng}}-W^{\mathrm{be}}|/W^{\mathrm{be}}$, where $W^{\mathrm{ng}}$ and $W^{\mathrm{be}}$ are the total ionization probabilities computed for the non-uniform grid and benchmark calculation, respectively.

Table 1 provides an overview of the calculation characteristics. As $h_{\max }$ increases, there is a corresponding decrease in the final number of grid points $N_g$ needed to describe the wave function, resulting in an enhanced computational efficiency $\eta$, which was estimated by comparing the propagation time between the non-uniform grid and benchmark calculations. Because we employ a dynamical grid, the efficiency $\eta$ does not scale as $1/N_g$. However, we observe a considerable improvement of about one order of magnitude even for this relatively short pulse. Additionally, we present the relative error on the total ionization probability, alongside with the maximum relative error on two energy intervals with different upper bonds for the energy: ${\mathcal{I}}_1=[0,0.3]$ and ${\mathcal{I}}_2=[0,0.5]$. As expected, the relative errors increase with larger $h_{\max }$ and higher electron energies.

We then repeated the calculation in VG using $dt=0.005$ and ${\ell }_{\max }=10$ to achieve converged results for the uniform grid. The obtained spectrum is depicted in Figure 2c, alongside spectra computed using non-uniform grids with various maximum mesh sizes. While the results are compared in the same gauge for consistency, however, we verified that the benchmark spectra in LG and VG are in excellent agreement with each other, with a maximum relative error on ${\mathcal{I}}_2$ of about $5\times {10}^{-3}$. A summary of the characteristics of the calculations in VG is presented in Table 2. It shows excellent accuracy and a significant gain in efficiency.

In the above calculations, we utilized an increment parameter $\alpha =0.05$. Experimentation with a smaller increment parameter, such as $\alpha =0.025$, revealed only a slight improvement in accuracy (not shown). Therefore, we generally recommend employing $\alpha =0.05$, as it is typically adequate to achieve the required level of accuracy in most practical applications.

The previous calculation, corresponding to a Keldysh parameter [27] $\gamma =2.13$, was relatively straightforward. However, it is also crucial to assess whether non-uniform grids can maintain high accuracy when electron tunneling and re-scattering events play significant roles. To address this aspect, we performed additional calculations using an 800 nm pulse with $N=6$ cycles at a peak intensity of ${10}^{14}$ W/cm², corresponding to $\gamma =1.06$. Converged results were achieved in VG for the uniform grid using $dt=0.0005$ and ${\ell }_{\max }=40$. In Figure 3a, we compare the spectra between the benchmark and the non-uniform grid calculations for $h_{\max }=0.25$. We observe excellent agreement in the reproduction of the complex spectrum structure across the entire energy interval $\mathcal{I}=[0,1]$, with the calculated maximum relative error ${ϵ}_r^{\max }\left(\mathcal{I}\right)=1.8\times {10}^{-2}$, while obtaining an efficiency gain of $\eta \approx 13$ using the non-uniform grid.

The gain in efficiency using the non-uniform grid increases towards the ratio $\eta \approx h_{\max }/h$ for longer pulse duration. To demonstrate this aspect, we performed additional calculations using a 150 nm pulse with $N=100$ cycles, at an intensity ${10}^{13}$ W/cm². The spectra, depicted in Figure 3b, again show excellent agreement in all features. In this case, the wave packet required the dynamical grid to extend to a final radius of 3500 a.u., leading to an efficiency gain by the non-uniform grid of $\eta \approx 20$.

We conclude this section by pointing out that the non-uniform grid method can be integrated with the t-SURFF surface-flux method [17]. In the t-SURFF approach, the energy-resolved ionization probability is translated into calculating the flux through a surface boundary positioned at a distance where the long-range Coulomb potential is assumed to be negligible. Typically, a surface boundary of approximately 200 atomic units is chosen, such that the non-uniform radial could only provide an additional improvement of a factor two to three. However, the non-uniform grid method offers several advantages on its own. Specifically, reproducing low-energy spectra in t-SURFF poses significant challenges and can be computationally demanding without sophisticated infinite-time propagation techniques [28,29], as the entire ionizing wave packet must cross the surface boundary. In contrast, the non-uniform grid has the ability to produce highly accurate spectra in the low-electron energy region. Furthermore, the t-SURFF method is not well-suited for computing HHG spectra or incorporating the dynamics of excited Rydberg states. In contrast, as will be demonstrated below, the non-uniform grid method naturally lends itself to handling such applications efficiently.

3.2. Application to High-Order Harmonic Generation

We now shift our focus to the computation of HHG spectra using the non-uniform grid, utilizing the same 800 nm pulse previously employed in Section 3.1. At each time step, we compute the time-dependent dipole moment $z\left(t\right)=\langle \mathsf{\Psi }\left(t\right)\left|z\right|\mathsf{\Psi }\left(t\right)\rangle$ using the modified Formula (27) for evaluating physical observables on the non-uniform grid. From the obtained value of $z\left(t\right)$, we compute the dipole acceleration $a\left(t\right)=d^{2}z\left(t\right)/d{t}^2$, which is then used to deduce the spectral density [30] as

$$S\left(\omega \right)=\frac{2}{3\pi c^3}{\left|\tilde{a}\left(\omega \right)\right|}^2,$$

(35)

where $\tilde{a}\left(\omega \right)=\int a\left(t\right)e^{-i\omega t}dt$ and c is the speed of light.

The benchmark calculation was performed in LG, using $dt=0.005$ and ${\ell }_{\max }=100$ to achieve converged results for a uniform grid with $h=0.01$. Subsequently, the calculation was repeated in LG for a non-uniform grid with $h_{\max }=0.5$ using the same time step and maximum angular momentum parameters. An additional calculation was performed in VG, using $dt=0.0005$ and ${\ell }_{\max }=40$.

The spectra resulting from these three calculations are shown in Figure 4, demonstrating remarkable agreement between the different gauges, despite the utilization of a maximum mesh size as large as $h_{\max }=0.5$. Detailed features across the entire spectra, including the structure near the plateau cutoff, are faithfully reproduced. The plateau cutoff is smooth for this relatively short pulse [31], and we also observe oscillations in the high-energy region due to the interference between recombining quantum trajectories. The excellent agreement can be attributed to the fact that the crucial part of the dipole moment resides near the atom, where the density of radial points remains high and the fact that returning trajectories do not acquire their highest kinetic energy during the far-range excursion of the electron. Nevertheless, for high-intensity, long-wavelength pulses, a reasonable maximum mesh size $h_{\max }$ needs to be chosen to achieve converged results. In a recent study [32], we used the non-uniform grid method to compute the HHG spectrum of Ne with a 2000 nm driving laser at an intensity of $2\times {10}^{14}$ W/cm² and obtained converged spectra with $h_{\max }=0.2$. This leads to a gain of efficiency of nearly $\eta =10$ when compared to the uniform grid. Consequently, the non-uniform grid method offers an accurate and efficient approach for calculating HHG spectra in atoms, including those generated by long-wavelength pulses with harmonics extending even into the water window [33]. The method will become particularly important for the treatment of high-Z targets (e.g., xenon).

3.3. Application to the Study of Rydberg States

The dynamics of highly excited Rydberg states in intense laser fields have sparked lively interest and debate. These states have been shown to survive high-intensity fields [3,4,34], can be used to create highly coherent wavepackets [35], and have also been suggested [36] to play a significant role in the low-energy electron peaks observed in the photoionization of atoms by mid-infrared pulses [37]. Moreover, Rydberg states with large angular momentum [38,39] represent attractive platforms for quantum information processing tasks due to their long lifetime and long-range dipole–dipole inteaction [40].

Describing highly excited Rydberg states presents challenges due to their large spatial extent. Computing Rydberg states up to principal quantum number $n=50$ necessitates a grid size exceeding 6000 atomic units. Utilizing a uniform grid with $h=0.01$, necessary for accurately capturing the rapid oscillations of the wave function near the atomic core, therefore involves diagonalizing the $N\times N$ field-free Hamiltonian (20) with $N=$ 60,000. This task is computationally demanding, even with state-of-the-art matrix diagonalization techniques. In contrast, employing a non-uniform radial grid with $h_{\max }=0.5$ only requires diagonalizing a matrix with $N=$ 15,000, a task that can be completed in about ten minutes on a desktop. Given that matrix diagonalization formally scales as $\mathcal{O}\left(N^3\right)$, a calculation using the uniform grid would take more than ten hours. While the calculation of Rydberg states up to $n=100$ could still be performed rather inexpensively using the non-uniform grid, it would currently be impractical on a uniform grid.

We computed Rydberg states ${\phi }_{n,\ell }$ for $\ell =0$ and $\ell =1$ up to $n=50$, and obtained a maximum relative error on their energy $E_{n,\ell }$ of less than $0.01\%$ for all calculated states. To illustrate the accuracy of the method, we employed a weak pulse with $N=20$ cycles at a peak intensity of ${10}^{12}$ W/cm², with two different central frequencies, $\omega =0.5$ and 0.52, near the ionization potential of hydrogen to both ionize the atom and populate highly excited Rydberg states via a one-photon transition. At the end of the pulse, we evaluated the population of Rydberg states $P_{n,\ell }={\left|\langle {\phi }_{n,\ell }|\mathsf{\Psi }\left(t\right)\rangle \right|}^2$ using the inner product rule (26). Following Li et al. [41], the energy-resolved ionization probability can be extended to the region of excited Rydberg states by defining a probability density $d{P}_n/dE=P_n(dn/dE)=P_n{n}^3$ for hydrogen since $E_n=-1/\left(2n^2\right)$. Using this prescription, we were able to obtain an accurate continuation of the ionization probability below the ionization threshold, as shown in Figure 5.

While our aim in this section was to demonstrate the accuracy and efficiency of the non-uniform grid and the validity of the inner-product formula to compute the population of excited states, such an approach opens the door to study strong-field processes with highly excited Rydberg states. This includes ionization suppression near channel closing [42], formation and dynamics of coherent Rydberg wave packets [35], interference stabilization effects [4,43], Freeman resonances [44,45], and low-energy electron peaks [37].

4. Conclusions and Outlook

We have presented a computational method to speed up frequently used single-active-electron computer codes for solving the time-dependent Schrödinger equation for the response of atoms to intense, short-pulse electromagnetic fields. The method was tested on three well-known benchmark problems involving the hydrogen atom, namely ejected-electron spectra in strong-field ionization, high-order harmonic generation, and the study of highly excited Rydberg states. Without losing any serious accuracy, speed-ups of at least an order of magnitude, but generally even more, were achieved. The Rydberg problem, in particular, would be intractable for states with principal quantum numbers significantly above $n=50$ with currently available computational resources if the usual approach with a uniform grid was employed. While other methods based on finite-element rather than finite-difference approaches are available, the non-uniform grid approach presented here is certainly competitive and should be attractive due to its relatively straightforward implementation.