Hydrogen Atom in Strong Elliptically Polarized Laser Fields within Discrete-Variable Representation

To study the dynamics of a hydrogen atom in an arbitrarily polarized laser field, we need to solve the corresponding TDSE (atomic units $m=e=\hbar=1$ are used throughout the text unless stated otherwise),

$\displaystyle i\frac{\partial}{\partial t}\psi({\bf r},t)=\left[H_{0}+V({\bf r},t)\right]\psi({\bf r},t)$

(1)

with the Hamiltonian $H_{0}=-\frac{1}{2}\Delta-\frac{1}{r}$ of the hydrogen atom in free space. We describe the interaction of the atomic electron with the laser field

$\displaystyle{\bf E}(t)=\frac{f(t)}{\sqrt{1+\epsilon^{2}}}E_{0}\left(\hat{\bf x}\cos\omega t+\hat{\bf y}\epsilon\sin\omega t\right)$

(2)

in the length gauge under the dipole approximation

$\displaystyle V({\bf r},t)=-{\bf r}\cdot{\bf E}(t)\,\,,$

(3)

where $E_{0}$, $0\leq\epsilon\leq 1$ and $\omega$ are the electric field amplitude, ellipticity and frequency of the laser field, respectively. The pulse envelope function is considered as

$\displaystyle f(t)=\cos^{2}(\frac{\pi t}{n_{T}T})~{}~{}~{},~{}~{}~{}-n_{T}T/2<t<n_{T}T/2\,\,,$

(4)

where the $n_{T}$ optical cycles of the period $T=2\pi/\omega$ are included during the pulse.

To solve TDSE, we follow the 2D DVR method [19, 20, 36, 37] where the electron wavefunction $\psi({\bf r},t)$ in the spherical coordinates $(r,\Omega)=(r,\theta,\phi)$ is expanded as

$\displaystyle\psi(r,\Omega,t)=\frac{1}{r}\sum_{j=1}^{N_{\Omega}}f_{j}(\Omega)u_{j}(r,t)$

(5)

over the basis $f_{j}(\Omega)$ associated with a 2D angular grid $\Omega_{j}$ ($j=1,2,...,N_{\Omega}$) on the unit sphere $\Omega$. In the works of V.S. Melezhik on 2D DVR in application to TDSE, the angular grids $\Omega_{j}=(\theta_{j_{\theta}},\phi_{j_{\phi}})$ were constructed as direct products of the nodes $\theta_{j_{\theta}}$ and $\phi_{j_{\phi}}$ of the 1D Gaussian quadratures over the $\theta$ and $\phi$ variables, respectively, where $j_{\theta}=1,...,N_{\theta}$, $j_{\phi}=1,...,N_{\phi}$ and $N_{\Omega}=N_{\theta}\times N_{\phi}$ [19, 20, 36]. This construction turned out to be very efficient in integration of the TDSE with three [19, 36, 38] and four [39, 40, 41] spatial variables. Particularly, with this approach the polarization of the high harmonics generated by hydrogen atoms were successfully investigated in elliptically polarized laser fields [19]. Recently [24], we have improved the computational scheme for the 3D stationary Schrödinger equation by implementing the 2D Popov and Lebedev nondirect product DVRs on the unit sphere $\Omega_{j}=(\theta_{j},\phi_{j})$ instead of the Gaussian 2D DVR. When integrating the stationary 3D Schrödinger equation for the hydrogen atom in crossed magnetic and electric fields, it was shown that the use of the 2D npDVR essentially accelerated the convergence of the computational scheme over $N_{\Omega}$ with respect to the 2D DVR scheme based on the Gaussian 1D quadratures with $N_{\Omega}=N_{\theta}\times N_{\phi}$ . This result motivated us to extend the 2D npDVR based on the Popov and Lebedev cubatures to the TDSE (1) for the hydrogen atom in a strong elliptically polarized laser field where convergence of the computational scheme becomes crucial due to complexity of the problem to be solved.

When angular grid points $\Omega_{j}=(\theta_{j_{\theta}},\phi_{j_{\phi}})$ are chosen as nodes of 1D Gaussian quadratures over the $\theta$ and $\phi$ variables, the basis functions $f_{j}(\Omega)$ are defined as

$\displaystyle f_{j}(\Omega)=\sum_{\nu=1}^{N_{\Omega}}\bar{Y}_{\nu}(\Omega)[Y^{-1}]_{\nu j}\,.$

(6)

Here the index $\nu$ numerates the pair $\{l,m\}$ ($m=-(N_{\phi}-1)/2,..,+(N_{\phi}-1)/2$ and $l=|m|,..,|m|+N_{\theta}-1$) related to the modified spherical harmonic

$\displaystyle\bar{Y}_{\nu}(\Omega)=\bar{Y}_{lm}(\Omega)=e^{im\phi}\sum_{l^{\prime}}{c_{l}^{l^{\prime}}}P_{l^{\prime}}^{m}(\theta)\,,$

(7)

where $P_{l^{\prime}}^{m}(\theta)$ are the associated Legendre polynomials. For the majority of $l$ and $l^{\prime}$ from a given above set, except the largest $l$ and $l^{\prime}$, $c_{l}^{l^{\prime}}=\delta_{ll^{\prime}}$ and $\bar{Y}_{lm}(\Omega)$ coincide with the spherical harmonics. The non zero $c_{l}^{l^{\prime}}$ for the largest $l$ and $l^{\prime}$ are found using the Gram-Schmidt orthogonalization procedure for $\bar{Y}_{lm}(\Omega)$ on the grid $\Omega_{j}$ [42]. The matrix $Y^{-1}$ is inverse to the $N_{\Omega}\times N_{\Omega}$ matrix $Y_{j\nu}=\sqrt{w_{j}}\bar{Y}_{\nu}(\Omega_{j})$, where $w_{j}=(2\pi{w_{j}}^{\prime})/N_{\phi}$ and ${w_{j}}^{\prime}$ are the weights of the Gaussian quadrature over the $\theta$ variable.

On the other hand, if the Popov or Lebedev cubature is used for generating the angular grid points $\Omega_{j}=(\theta_{j},\phi_{j})$ with $j=1,..,N_{\Omega}$, then $f_{j}(\Omega)$ is defined as

$\displaystyle f_{j}(\Omega)=\sum_{\nu=1}^{N_{\Omega}}\sqrt{4\pi w_{j}}\Phi_{\nu}(\Omega)\Phi^{*}_{\nu}(\Omega_{j})\,\,,$

(8)

where $w_{j}$ is the Popov or Lebedev weight ($\sum_{j=1}^{N_{\Omega}}{w_{j}}=1.0$) (for Lebedev and Popov grid points and weights see [43] and [25, 26, 27] , respectively). The function $\Phi_{\nu}$ is a special linear combination of the spherical harmonics, $\Phi_{\nu}(\Omega)=\sum_{\mu}S_{\nu\mu}Y_{\mu}(\Omega)$ obtained within the method described in our previous work [24]. Here $\hat{S}$ refers to a $N_{\Omega}\times\bar{N}_{\Omega}$ matrix ( with $\mu=1,...,\bar{N}_{\Omega}\geq N_{\Omega}$ ) which makes the set of $N_{\Omega}$ basis functions $\Phi_{\nu}(\Omega)$ orthonormal over the grid points $\Omega_{j}$ in the sense of the Popov or Lebedev cubatures.

Thus, the problem is reduced to a system of Schrödinger-type equations which should be solved for the unknown $N_{\Omega}$-dimensional vector ${\bf u}(r,t)=\{\sqrt{w_{j}}u_{j}(r,t)\}_{1}^{N_{\Omega}}=\{\sqrt{w_{j}}r\psi(r,\Omega_{j},t)\}_{1}^{N_{\Omega}}$,

$\displaystyle i\frac{\partial}{\partial t}{\bf u}(r,t)=\left[\hat{H}_{0}(r)+\hat{V}(r,t)\right]{\bf u}(r,t)\,\,,$

(9)

where the elements of the $N_{\Omega}\times N_{\Omega}$ matrices $\hat{H}_{0}(r)$ and $\hat{V}(r,t)$ are represented as

	$\displaystyle H_{0jj^{\prime}}(r)$	$\displaystyle=$	$\displaystyle-\frac{1}{2}(\delta_{jj^{\prime}}\frac{d^{2}}{dr^{2}}-\frac{L_{jj^{\prime}}^{2}}{r^{2}})-\frac{1}{r}\delta_{jj^{\prime}}\,,$		(10)
	$\displaystyle V_{jj^{\prime}}(r,t)$	$\displaystyle=$	$\displaystyle V(r,\Omega_{j},t)\delta_{jj^{\prime}}\,.$		(11)

In the Gaussian case, the elements of the ${\hat{L}}^{2}$ matrix are defined as

$\displaystyle L_{jj^{\prime}}^{2}=\sum_{\nu=1}^{N_{\Omega}}{\sqrt{w_{j}w_{j^{\prime}}}}\bar{Y}_{\nu}(\Omega_{j})l(l+1){\bar{Y}_{\nu}}^{\ast}(\Omega_{j^{\prime}})\,.$

(12)

whereas in the Popov and Lebedev cases, they are written as

$\displaystyle L_{jj^{\prime}}^{2}=4\pi\sum_{\nu=1}^{N_{\Omega}}{\sqrt{w_{j}w_{j^{\prime}}}}\sum_{\mu=1}^{\bar{N}_{\Omega}}S_{\nu\mu}l_{\mu}(l_{\mu}+1){{Y}_{\mu}}(\Omega_{j}){\Phi}_{\nu}^{\ast}(\Omega_{j^{\prime}})\,,$

(13)

with $\bar{N}_{\Omega}\geq N_{\Omega}$ and $\mu=\{l_{\mu},m_{\mu}\}$ defined in [24].

To propagate the wavefunction in time $t_{n}\rightarrow t_{n+1}=t_{n}+\Delta t$, the component-by-component split-operator method [44] is applied as

$\displaystyle{\bf u}(r,t_{n}+\Delta t)=\exp[-\frac{i}{2}\Delta t\hat{V}(r,t_{n})]\exp[-i\Delta t\hat{H}_{0}(r)]\exp[-\frac{i}{2}\Delta t\hat{V}(r,t_{n})]{\bf u}(r,t_{n})\,\,.$

(14)

Since the interaction potential $\hat{V}(r,t)$ is diagonal in our 2D DVRs over both angular grid distributions (either Gaussian $\Omega_{j}=(\theta_{j_{\theta}},\phi_{j_{\phi}})$ or Popov (Lebedev) $\Omega_{j}=(\theta_{j},\phi_{j})$), the first and last steps of the procedure(14) represent simple multiplications of the vectors $\exp\left[-(i/2)\Delta tV(r,\Omega_{j},t)\right]$ and ${\bf u}(r,t)$ which results in ${\bf u}(r,t+\frac{1}{4}\Delta t)$. The intermediate step is performed after diagonalizing the nondiagonal part (${\hat{L}}^{2}$ operator) which transforms $\hat{H}_{0}$ to a diagonal matrix $\hat{A}$. Then we can approximate its related exponential expression according to

$\displaystyle\exp[-i\Delta t\hat{A}]\approx{\left(1+\frac{i}{2}\Delta t\hat{A}\right)}^{-1}\left(1-\frac{i}{2}\Delta t\hat{A}\right)$

(15)

which ensures the desired accuracy of the numerical algorithm (14).

Thus, we reach the following initial value problem for the intermediate step:

$\displaystyle\left(1+\frac{i}{2}\Delta t\hat{A}\right)\bar{\bf u}(t_{n}+\frac{3}{4}\Delta t)=\left(1-\frac{i}{2}\Delta t\hat{A}\right)\bar{\bf u}(t_{n}+\frac{1}{4}\Delta t)\,,$

(16)

$\displaystyle\bar{\bf u}(r=0,t_{n}+\frac{3}{4}\Delta t)=\bar{\bf u}(r=r_{m},t_{n}+\frac{3}{4}\Delta t)=0~{}~{},~{}~{}r_{m}\rightarrow\infty\,\,.$

(17)

Here the transformed wavefunction $\bar{\bf u}$ and operator $\hat{A}$ are defined based on the transformation method which is used for diagonalizing the ${\hat{L}}^{2}$ operator. When using Gaussian grid points, we apply a simple unitary transformation $M_{j\nu}=\sqrt{w_{j}}\bar{Y}_{\nu}(\Omega_{j})$ [20, 36], then

$\displaystyle\bar{\bf u}=\hat{M}{\bf u}\,,$

(18)

$\displaystyle\hat{A}_{\nu\nu^{\prime}}=\left(\hat{M}\hat{H}_{0}\hat{M^{\dagger}}\right)_{\nu\nu^{\prime}}=\left[-\frac{1}{2}\frac{d^{2}}{dr^{2}}+\frac{l(l+1)}{2r^{2}}-\frac{1}{r}\right]\delta_{\nu\nu^{\prime}}$

(19)

with $\nu=\{l,m\}$.

In the Popov and Lebedev cases, we use the conventional numerical diagonalization method. We find the eigenvectors and eigenvalues of the ${\hat{L}}^{2}$ matrix in these representations. Then the transformation matrix $\hat{W}$ is constructed from the numerically calculated eigenvectors which are stored on its columns so that $\hat{W}^{-1}{\hat{L}}^{2}\hat{W}=\hat{J}$, where $\hat{J}$ is the diagonal matrix of the acquired eigenvalues $J_{\nu}$. Thus $\bar{\bf u}$ and $\hat{A}$ are represented as

$\displaystyle\bar{\bf u}=\hat{W}^{-1}{\bf u}\,,$

(20)

$\displaystyle\hat{A}_{\nu\nu^{\prime}}=\left(\hat{W}^{-1}\hat{H}_{0}\hat{W}\right)_{\nu\nu^{\prime}}=\left[-\frac{1}{2}\frac{d^{2}}{dr^{2}}+\frac{J_{\nu}}{2r^{2}}-\frac{1}{r}\right]\delta_{\nu\nu^{\prime}}$

(21)

with $\nu=1,2,...,N_{\Omega}$.

Finally, before applying the last operation in (14), $\bar{\bf u}$ should be transformed to ${\bf u}$ through the operation ${\bf u}=\hat{M^{\dagger}}\bar{\bf u}$ (${\bf u}=\hat{W}\bar{\bf u}$) while using Gaussian (Popov or Lebedev) quadratures (cubatures).

It should be emphasized that both transformation matrices (either $\hat{M}$ or $\hat{W}$) are independent of $r$ and $t$.

For approximating the Schrödinger equation (9) over the radial variable $r$ the six-order finite-difference approximation is used on a quasi-uniform grid $\{r_{j}\}_{1}^{N_{r}}$ on the interval $r\in[0,r_{m}]$ constructed by mapping the $r$ variable to $x\in[0,1]$ with the formula $r=r_{m}[\exp(4x)-1]/[\exp(4)-1]$. To avoid the artificial reflection of the electron at the radial grid boundary, the wave function is multiplied after each time step by the absorbing “mask function” suggested in Ref.[45]. The total number of operations at each time step (14) in our algorithm is $N_{r}(2N_{\Omega}^{2}+\alpha N_{\Omega})$ where $\alpha=(2\times 7)^{2}\simeq 200$ is given by the width of the band of the band matrices $\hat{A}$ in (19) and (21) with the finite-difference approximation of the radial derivatives with respect to the variable $r$. Here we have an obvious advantage with respect the Tong computational scheme [16] for numerical integration of the Shrödinger equation (1) where the number of operations at each time step is $N_{r}N_{\theta}N_{\phi}(N_{r}+N_{\theta}+N_{\phi})$, i.e. quadratically dependent on the number of radial grid points $N_{r}$. The computational time of our scheme linearly depends on $N_{r}$ and $\sim N_{\Omega}^{\gamma}$ with $1\lesssim\gamma<2$ approaching 1 for not very big $N_{\Omega}$ [20]. Linear dependence of the computational time on $N_{r}$ permits computations on the radial grids up to $r_{m}\ =1000$ and $\Delta x=1/N_{r}=0.00025$. The time step of integration is chosen to be $\Delta t=0.05$ ($1/2200$ of one optical cycle).

Once the wavefunction $\psi({\bf r},t)$ is calculated, one can investigate the physical properties of the atom in the laser field such as photoelectron momentum spectra, photoexcitation or HHG etc. Here we focus on the ionization and excitation probabilities.

We calculate the transition probabilities in the hydrogen atom due to interaction with a laser pulse by projecting the electron wave-packet at the end of the pulse $\psi({\bf r},t=n_{T}T)$ on the bound

$\displaystyle P_{nlm}=|\langle\phi_{nlm}|\psi\rangle|^{2}=|d_{nlm}|^{2}$

(22)

and continuum

$\displaystyle\frac{dP_{lm}}{dE}(E)=|\langle\phi_{klm}|\psi\rangle|^{2}=|c_{klm}|^{2}$

(23)

states of the atomic spectrum in free space. The hydrogen eigenfunctions in the above formulas (22) and (23) of the bound

$\displaystyle\phi_{nlm}({\bf r})=R_{nl}(r)Y_{lm}(\Omega)$

(24)

and continuum

$\displaystyle\phi_{klm}({\bf r})=\sqrt{\frac{2}{\pi k}}\frac{F_{l}(\eta,kr)}{r}Y_{lm}(\Omega)$

(25)

spectrum form an orthonormal and complete basis set which can be used for expansion of the electron wave-packet $\psi({\bf r},t)$ after interaction with the laser pulse

$\displaystyle\psi({\bf r},t=n_{T}T)=\sum_{nlm}d_{nlm}\phi_{nlm}({\bf r})+\sum_{lm}\int_{0}^{+\infty}c_{klm}\phi_{klm}({\bf r})dE\,\,,$

(26)

where $E$ denotes the photoelectron energy, $k=\sqrt{2E}$, $\eta=-1/k$ and $F_{l}(\eta,kr)$ is the Coulomb wavefunction regular as $kr\rightarrow 0$ .

Thus, the population of the $n$-th bound state and the ATI spectrum of the hydrogen atom after interaction with the laser pulse are described as

$\displaystyle P_{n}=\sum_{l=0}^{n-1}\sum_{m=-l}^{l}P_{nlm}$

(27)

and

$\displaystyle\frac{dP}{dE}(E)=\sum_{l=0}^{l_{m}\rightarrow\infty}\sum_{m=-l}^{l}\frac{dP_{lm}}{dE}(E)\,.$

(28)

Together with the above equations we also use the relation between the population probabilities

$\displaystyle\sum_{n=1}^{\infty}P_{n}+\int_{0}^{+\infty}\frac{dP}{dE}dE=\langle\psi|\psi\rangle\approx 1\,.$

(29)

Due to the action of the absorbing “mask function” at the edge of radial integration $r_{m}$, in our computation we have to take into account the effect that the normalization integral $\langle\psi|\psi\rangle$ decreases in time and deviates from unit. Then the total ionization yield is obtained by

$\displaystyle P_{ion}=\int_{0}^{+\infty}\frac{dP}{dE}dE+\left(1-\langle\psi|\psi\rangle\right)\,,$

(30)

and $\sum_{n=1}^{\infty}P_{n}=P_{g}+P_{ex}$, where $P_{g}$ and $P_{ex}$ denote the ground state population and total excitation probability.

To circumvent the seemingly intractable problem of calculating the populations $P_{n}$ as $n\rightarrow+\infty$ , we propose an alternative method which does not require numerical analysis for very large $n,l,m$ values. It is based on the known property of the Coulomb wave functions $\phi_{nlm}({\bf r})$ (24) and $\phi_{klm}({\bf r})$ (25): the functions $n^{3/2}\phi_{nlm}({\bf r})$ and $\phi_{klm}({\bf r})$ transform one to another as $n\rightarrow+\infty$ and $k\rightarrow 0$ with the replacement $n\leftrightarrow 1/(ik)$ [46, 47].

In this method, we stop calculating $P_{n}$ at some rather large $n=N^{\prime}$

$\displaystyle P_{ex}=\sum_{n=2}^{\infty}P_{n}=\sum_{n=2}^{N^{\prime}}P_{n}+\sum_{n=N^{\prime}+1}^{\infty}P_{n}$

(31)

and evaluate the remaining sum $\sum_{n=N^{\prime}+1}^{\infty}P_{n}$ in the above summation according to

$\displaystyle\sum_{n=N^{\prime}+1}^{+\infty}P_{n}\simeq\int_{n=N^{\prime}+1}^{+\infty}P_{n}dn=\int_{N^{\prime}+1}^{+\infty}n^{3}P_{n}d(-\frac{1}{2n^{2}})=\int_{E^{\prime}}^{0}\frac{dP}{dE}(E<0)dE\,,$

(32)

where $E^{\prime}=-\frac{1}{2(N^{\prime}+1)^{2}}$ , and $n^{3}P_{n}=\frac{dP}{dE}(E<0)$. It is well known [48, 49] (see also [17]) that in the limit $n\rightarrow+\infty$ the oscillator strength density $\frac{dP}{dE}(E<0)$ approaches the value $\frac{dP}{dE}(-0)$ coinciding with the limit $\frac{dP}{dE}(+0)$ of the population of the states of the continuum spectrum at $k\rightarrow 0(E\rightarrow+0)$.

Thus, by calculating $\frac{dP}{dE}(E>0)$ from Eq.(28) at some positive small values of $E~{}(E\rightarrow+0)$, where small $l,m$ give a principal contribution, we can find a trend curve $f(E)=\frac{dP}{dE}(0)+AE+BE^{2}+...$ passing through these points and the points $n^{3}P_{n}=\frac{dP}{dE}(E<0)$ calculated on the negative side of $E$ which belong to the $n=N^{\prime},~{}N^{\prime}-1,...$ states. With the calculated $f(E)$ we perform infinite summation in (32)

$\displaystyle\sum_{n=N^{\prime}+1}^{+\infty}P_{n}=\sum_{n=N^{\prime}+1}^{n_{m}}\frac{f(E_{n})}{n^{3}}+\int_{E_{1}}^{0}f(E)dE$

(33)

where $E_{1}=E_{n_{m}+1}$ and $n_{m}$ is chosen to be as large as $30$ after which the difference between summation over $n$ and integration over $E$ becomes negligible. Finally, we can obtain $P_{ex}$ by Eq.(31) and then the total ionization yield by $P_{ion}=1-P_{g}-P_{ex}$.

First, we consider a hydrogen atom in a laser field of linear polarization ($\epsilon=0$). This case permits the separation of the azimuthal angle $\phi$ if the electric field is oriented along the $z$-axis and the 2D DVR basis constructed with the Gaussian scheme (6) reduces to the 1D DVR over the $\theta$ variable with $N_{\Omega}=N_{\theta}$ ($l=0,1,..,N_{\Omega}-1$, and $m=0$). Such simplification of the problem gives us the possibility to get the most accurate result and use it as a test for the convergence of the 2D npDVRs constructed with the Popov and Lebedev cubatures where the separation of angular variables are not allowed, regardless of the chosen orientation for the electric field in space. In Table 1, we demonstrate the convergence of the scheme based on the Gaussian 1D DVR with respect to the number of the basis function $N_{\Omega}$ in calculating the populations $P_{n}$ of the ground ($n=1$) and excited states up to $n=10$. It is shown that including $N_{\Omega}=35$ basis functions in the 1D DVR Gaussian scheme (where $l\leq 34$ and $m=0$) fixes four significant digits after the decimal point in $P_{1}$ and gives the relative accuracy $\sim 9\%$ in calculating $P_{n}$ for $2\leq n\leq 10$. When increasing the number of basis functions to $N_{\Omega}=95$ ($l\leq 94$, $m=0$), eight significant digits are fixed in $P_{1}$, and the relative accuracy reaches the value $\sim 0.5\%$ in $P_{n}$ for $2\leq n\leq 10$. The distribution of the populations $P_{nl}$ over the angular momentum $l$ is presented in Fig.1. Here, the most probable angular momentum states are observed with $l\leqslant 5$.

In Fig.2, we demonstrate the convergence $N_{\Omega}\rightarrow+\infty$ of the transition probabilities $\frac{dP}{dE}(E)$ to the states of the continuum spectra of the hydrogen atom due to interaction with a laser pulse linearly polarized along the $z$-axis. The calculations are performed with the 1D DVR Gaussian scheme where the $\phi$-variable is separated and $m=0$. In the calculations, the highest angular momentum in Eq.(28) was fixed at $l_{m}=35$ by the requirement that the error in cutting off summation (28) does not exceed the error in using the limited DVR basis. The spectrum contains several peaks whose position does not depend on the number of basis functions $N_{\Omega}$ and are almost placed at the peak positions $E_{p}$ of the energy spectrum in a multiphoton process [50, 51, 52]

$\displaystyle E_{p}=p\hbar\omega-U_{pe}-I_{i}\,\,,$

(34)

where $I_{i}=13.6$ eV is the ionization energy of the hydrogen atom, $\hbar\omega=1.55$ eV is the photon energy, and $U_{pe}=I/\left(4\omega^{2}\right)=5.967$ eV is the ponderomotive energy of our laser field with intensity $I=10^{14}$ W/cm². The calculated $\frac{dP}{dE}(E)$ is qualitatively in agreement with Ref.[51] where a $25$ fs laser pulse was used. For this high intensity, the first expected peak corresponding to the 13-photon ionization, $p=13$ in Eq.(34), is splitted, just as in [51]. The various mechanisms leading to the doubling of the first ATI peak was intensively discussed earlier [33, 53, 54, 55]. In the considered case, apparently, the intermediate resonant population of some of the levels dressed by strong laser field is responsible for this effect [33]. However, to clarify this issue, a detailed calculation of the level deformations by a strong laser pulse is required.

The distribution $\frac{dP_{l}}{dE}(E)$ over electron angular momentum $l$ in the continuum for different energies $E$ between successive peaks ($E=E_{p}+0.01$ (in a.u.), where $p$ corresponds to $13,14,15,16-$photon ionization) is presented in Fig.3. As can be seen in this figure, by increasing the electron energy, the contribution of higher $l$-values to the continuous spectrum is increased. However, the main contribution comes from angular momenta $0\leqslant l\leqslant 12$.

In Table 2, the total ionization yield $P_{ion}$ which is calculated by integrating the transition probabilities $\frac{dP}{dE}(E)$ over all possible states of the ionized electron in its continuum (30), and the total excitation probability $P_{ex}=1-P_{g}-P_{ion}$ are presented. These results demonstrate the convergence of the computational scheme in calculating $P_{ion}$ and $P_{ex}$ at the accuracy level of the order of the fourth significant digit after the decimal point.

The total excitation probability $P_{ex}$ can also be calculated by direct summation of $P_{nl}$ over $n$ and $l$. However, this requires taking into account the contribution of large $n$ and $l$ values which in some cases can be considerable [34, 56, 57] and greatly complicates the problem computationally. Therefore, we propose an alternative method described in the previous section , equations (31) and (33), to calculate $P_{ex}$ without summation over large $n$ and $l$. To get the probability $P_{ex}$ according to these formulas (the procedure is illustrated in Fig.4), we have calculated the populations $P_{n}$ of the hydrogen atom bound states up to $n=7$ , according to (27), and followed the method described in Sec.2.2 for calculating the populations of states with $n\geq 8$. The points $\frac{dP}{dE}(E_{n})=n^{3}P_{n}$ on the negative side of the electron energy in Fig.4 belong to $n=6$ and $n=7$ states. On the positive side, three $E$ points were chosen very close to the zero energy where $\frac{dP}{dE}(E)$ were calculated by Eqs.(23,28). By these five points the polynomial function $f(E)$ was fixed for calculating $\sum_{n=8}^{+\infty}P_{n}$ , according to the relation (33). In Fig.4 we see a good agreement between the calculated curve $f(E)$ and the corresponding values calculated for $n$ up to $n=20$ by the direct way (22). It confirms that the choice of $n=6$ and $n=7$ when constructing $f(E)$ gives a good approximation which can be used to evaluate $P_{n}$ for large $n$. The probabilities $P_{ex}$ calculated with this method and $P_{ion}=1-P_{g}-P_{ex}$ as a function of $N_{\Omega}$ are presented in Table 3. By comparing the results presented in Tables 2 and 3, one can conclude that our method (31,33) proposed for calculating the total excitation probability $P_{ex}$ and ionization yield $P_{ion}$ generates values which coincide, within a given accuracy, with the corresponding values obtained above by integrating the electron transition probabilities over the states of an electron in the continuum spectrum.

The values of $P_{g}$, $P_{ex}$ and $P_{ion}$, calculated above with high accuracy as $N_{\Omega}=N_{\theta}\rightarrow+\infty$ in the framework of the Gaussian 1D DVR, hereinafter are used as “exact” reference results to investigate the convergence and evaluate the accuracy of the computational schemes based on the Popov and Lebedev 2D npDVRs, where the separation of the $\phi$-variable is not allowed regardless of the laser field orientation. In Fig.5 , we demonstrate the convergence of the computational schemes based on 2D npDVRs by the example of calculating the population of the ground state $P_{g}$ as a function of the number of the angular grid points on the unit sphere $N_{\Omega}$ (the number of basis functions (8) in 2D npDVR expansion (5)). As can be seen in this figure, the ground state population, calculated with the Popov or Lebedev schemes, tends to the “exact” value obtained according to the Gaussian scheme (yellow dashed line). We have also performed the calculation of $P_{g}(N_{\Omega})$ for the nonseparable case of the Gaussian 2D DVR, when the laser field is oriented along the $x$-axis and $N_{\Omega}=N_{\theta}\times N_{\phi}$. Figure 5 displays much faster convergence of the Popov and Lebedev schemes compared to the Gaussian 2D DVR. The slowdown of the latter can be explained by the known effect of clustering of grid points $\Omega_{j}$ in the Gaussian scheme near the poles $z=0$ as $N_{\Omega}\rightarrow\infty$ [58] and a more adequate approximation of the angular part of the 3D TDSE in 2D npDVRs.

A comparison of the different applied schemes in calculating $P_{g}$, $P_{ex}$ and $P_{ion}$ is summarized in Table 4. Here the largest number of employed Popov and Lebedev grid points is $N_{\Omega}=672$ and $1202$, respectively. A comparison of the most accurate values obtained with the Lebedev 2D npDVR scheme at $N_{\Omega}=1202$ with the “exact” values calculated using the 1D Gaussian DVR for field polarization along the $z$-axis demonstrates the achieved absolute accuracy of this scheme of the order of $\sim 10^{-5}$ for $P_{g}$ and $\sim 10^{-4}$ for $P_{ex}$ and $P_{ion}$. This estimated accuracy also coincides with the value which is obtained through the generally accepted procedure for calculating the error of the numerical methods on a sequence of condensing grids where the error is equal to the attained minimum deviation of the calculated value on two nearest most detailed grids. Indeed, in our case, for the deviation of the populations calculated on the most detailed grid from its closest number of nodes we get $\Delta P_{g}=\mid P_{g}(N_{\Omega}=974)-P_{g}(N_{\Omega}=1202)\mid\simeq 2\times 10^{-5}$, $\Delta P_{ex},\Delta P_{ion}\simeq 10^{-4}$. This method is also used in the next section for evaluating the accuracy of calculations. The above results were obtained on a radial grid with $r_{m}=500$ and $\Delta x=0.0005$. To evaluate the error of the approximation over radial variable we have performed calculations with increasing radial box size and decreasing step of integration up to $r_{m}=1000$ and $\Delta x=1/N_{r}=0.00025$. The results given in Table 5 demonstrate that the approximation error over the radial variable on the grid with $r_{m}=500$ and $\Delta x=0.0005$ is less than due to truncation of the 2D npDVR. Therefore, hereafter we choose this radial grid. The step of integration $\Delta t=0.05$ ($1/2200$ of one optical cycle) over time, which we use, also gives an error that is significantly less than the error introduced by truncating the 2D npDVR.

Although our results deviate from the result of Ref.[17] (the maximal deviation is for $P_{ex}$ values which differ from one another twice), our three DVR computational schemes converge to the same value in the frame of the reached accuracy. Moreover, the deviation of our results from the values obtained in [17] exceeds our accuracy by an order of magnitude.

In the next section we apply the 2D npDVRs for a hydrogen atom in elliptically polarized laser fields in the entire range of possible change in ellipticity.

The dynamics of a hydrogen atom in a laser field of three different polarizations $\epsilon=0.25$, $0.5$, $1.0$ has been investigated with the developed Popov and Lebedev 2D npDVRs. As was mentioned above, the 1D Gaussian DVR is not applicable for a hydrogen atom in a laser field with $\epsilon\neq 0$. The calculated ground state probabilities $P_{g}$ are presented in Fig.6, which demonstrates the convergence of our 2D DVR computational schemes as $N_{\Omega}\rightarrow+\infty$. As in the case of linear polarization, the Popov and Lebedev 2D npDVRs schemes demonstrate more rapid convergence compared to the 2D Gaussian DVR for circular polarization as well. As in the case of linear polarization, all three computational schemes converge to the same values, which, however, significantly deviate from the results obtained in [17]. This deviation decreases with increasing ellipticity and becomes minimal at $\epsilon=1$.

The distributions of the population ($P_{n}$ and $P_{nl}$) of excited states over principal quantum number $n$ and angular momentum $l$ are presented in Fig.7. While for ellipticities $\epsilon=0.0,0.25,0.5$ a maximum is observed at $n=5$ bound state and a strong decrease in $P_{n}$ for higher excited states (in agreement with [34]), the $n$ distribution shows a different trend for the circularly polarized laser field.

To calculate the probabilities of total excitation $P_{ex}$ and ionization $P_{ion}$ of the hydrogen atom, depending on the ellipticity of the laser field, we have used the Lebedev 2D npDVR scheme, which allowed us to carry out calculations with the largest number of basis DVR functions up to $N_{\Omega}=1202$. In this case, to approximate the infinite summation (32) when calculating the probability $P_{ex}$ of total excitation (31), we use the computational procedure (31), (32), whose efficiency was demonstrated in Sec.3.1 for linear polarization of the laser field. Figure 8 illustrates this computational procedure for $\epsilon=0.5$ and $1$. The ground state population $P_{g}$ and thus calculated probabilities $P_{ex}$ and $P_{ion}$ are presented in Tables 6, 7 and 8. Figures 9 also depict the obtained $P_{ex}$ and $P_{ion}$ values. The data presented in Tables 6, 7 and 8 and in Figs.9 convincingly demonstrate the convergence of the Lebedev 2D npDVR scheme as $N_{\Omega}\rightarrow+\infty$ and give the estimate $\Delta P_{x}\simeq\mid P_{x}(N_{\Omega}=974)-P_{x}(N_{\Omega}=1202)\mid$ of the achieved absolute accuracy in calculating $P_{g}$, $P_{ex}$ and $P_{ion}$, which is improved from the values of the order $\sim 10^{-4}$ at $\epsilon=0$ to $10^{-6}$ at $\epsilon=1$. It should be noted that the errors due to the finiteness of the radial box size $r_{m}$ and the steps of integration over the radial and time variables do not exceed this achieved accuracy. Table 9 illustrates that the radial grid with $r_{m}=500$ and $\Delta x=1/N_{r}=0.0005$ already meets this requirement when calculating $P_{ex}$ and $P_{ion}$. The difference between the calculated values $P_{g}$, $P_{ex}$ and $P_{ion}$ and most accurate values [17] calculated so far with alternative methods decreases with increasing ellipticity to $\epsilon=1$. However, in the entire range of $\epsilon$ variation, this difference remains an order of magnitude greater than the achieved accuracy of our calculations. We explain this deviation by the fact that we managed to significantly improve, compared to [17], the approximation of the angular part of the 3D TDSE using our 2D npDVR, which ensured fast convergence, and high computational efficiency of our algorithm where the computational time is proportional to $\sim N_{r}(2N_{\Omega}^{2}+200N_{\Omega})$. Perhaps, this is also due to the fact that in the pseudospectral method of [17] instead of the unperturbed bound states of the hydrogen atom $R_{nl}(r)Y_{lm}(\theta,\phi)$ , pseudostates are used starting from $n>16$, although there is a fairly large set of them.