Photoionization of few electron systems: a hybrid coupled channels approach

We discretize our N-electron wavefunction by channel wave functions chosen as anti-symmetrized products of ionic states with a numerical one-electron basis. To this we add the wave function of the neutral ground state, resulting in the expansion

$$\begin{eqnarray}&&|\Psi (t)\rangle \approx \mathop{\sum }\limits_{\mathcal{I}}|\mathcal{I}\rangle {{C}_{\mathcal{I}}}(t)+|\mathcal{G}\rangle {{C}_{\mathcal{G}}}(t),\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&|\mathcal{I}\rangle =\mathcal{A}[|i\rangle |I\rangle ].\end{eqnarray} \tag{ 6 }$$

Here, $\mathcal{A}$ indicates anti-symmetrization, $|i\rangle$ are functions from a numerical one-electron basis, $|I\rangle$ and $|\mathcal{G}\rangle$ are (N-1) and N particle functions respectively and ${{C}_{\mathcal{I}}}(t)$, ${{C}_{\mathcal{G}}}(t)$ are the time dependent coefficients.

For the single electron basis $|i\rangle$ we use a finite element representation on the radial coordinate times spherical harmonics on the angular coordinates

$$\begin{eqnarray}&&|i(\vec{r})\rangle =|{{f}_{i}}(r)\rangle |{{Y}_{{{l}_{i}}{{m}_{i}}}}(\hat{r})\rangle .\end{eqnarray} \tag{ 7 }$$

On each finite element we use high order scaled Legendre polynomials as basis functions. The typical orders we use are 10–14. The details of the finite element approach used here can be found in [23, 24]. A brief description is given in appendix for the convenience of the reader. We refer to this basis as the active electron basis.

We choose $|I\rangle$ to be the eigenstates of the single ionic Hamiltonian obtained from the multi-reference configuration interaction singles doubles (MR-CISD) level of quantum chemistry. $|\mathcal{G}\rangle$ is chosen as the ground state of the system, also obtained from the MR-CISD level of quantum chemistry. These quantum chemistry wavefunctions are constructed with an atom centered primitive Gaussian basis as the starting point. While $|\mathcal{G}\rangle$ is the lowest eigenvector of the N particle Hamiltonian as obtained from COLUMBUS, it is not the ground state of the Hamiltonian in our basis: by treating one of the electrons with the active electron basis that is superior to the Gaussian basis one further improves the ground state energy.

The wavefunctions $|\mathcal{G}\rangle$ and $|I\rangle$ can be represented in a general form as sums of determinants:

$$\begin{eqnarray}&&|I\rangle =\mathop{\sum }\limits_{{{p}_{1}},{{p}_{2}},..,{{p}_{n-1}}}{{d}_{{{p}_{1}},{{p}_{2}},..,{{p}_{n-1}}}}|\mathcal{A}[{{\phi }_{{{p}_{1}}}}{{\phi }_{{{p}_{2}}}}...{{\phi }_{{{p}_{n-1}}}}]\rangle \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&|\mathcal{G}\rangle =\mathop{\sum }\limits_{{{p}_{1}},{{p}_{2}},..,{{p}_{n}}}{{d}_{{{p}_{1}},{{p}_{2}},..,{{p}_{n}}}}|\mathcal{A}[{{\phi }_{{{p}_{1}}}}{{\phi }_{{{p}_{2}}}}...{{\phi }_{{{p}_{n}}}}]\rangle ,\end{eqnarray} \tag{ 9 }$$

where ${{\phi }_{k}}$ are the Hartree–Fock orbitals of the neutral system. The same set of Hartree–Fock orbitals are used to construct both ionic and neutral CI functions. This allows us to use simple Slater–Condon rules to compute any matrix elements between them.

The explicit inclusion of the neutral ground state is motivated by the fact that, while the ionization process itself may be well described by one or a few ionic channels, the initial ground state may be more strongly correlated. In order to avoid inclusion of many ionic states just to describe the initial state, we include the neutral ground state explicitly, thereby reducing the number of basis functions needed. This idea can be easily extended to include any specific correlated state that is of importance to a particular process.

Substituting the ansatz (5) into the TDSE (1) yields a set of coupled ordinary differential equations for the time dependent coefficients:

$$\begin{eqnarray}&&i\left[ \langle \mathcal{G}|\mathcal{G}\rangle \frac{{\rm d}{{C}_{\mathcal{G}}}}{{\rm d}t}+\langle \mathcal{G}|\mathcal{I}\rangle \frac{{\rm d}{{C}_{\mathcal{I}}}}{{\rm d}t} \right]=\langle \mathcal{G}|\hat{H}|\mathcal{G}\rangle {{C}_{\mathcal{G}}}+\langle \mathcal{G}|\hat{H}|\mathcal{I}\rangle {{C}_{\mathcal{I}}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&i\left[ \langle \mathcal{I}|\mathcal{G}\rangle \frac{{\rm d}{{C}_{\mathcal{G}}}}{{\rm d}t}+\langle \mathcal{I}|\mathcal{I}\rangle \frac{{\rm d}{{C}_{\mathcal{I}}}}{{\rm d}t} \right]=\langle \mathcal{I}|\hat{H}|\mathcal{G}\rangle {{C}_{\mathcal{G}}}+\langle \mathcal{I}|\hat{H}|\mathcal{I}\rangle {{C}_{\mathcal{I}}}.\end{eqnarray} \tag{ 11 }$$

The time-derivative of the coefficient vector is multiplied by an overlap matrix composed of the blocks $\langle \mathcal{G}|\mathcal{G}\rangle$, $\langle \mathcal{G}|\mathcal{I}\rangle$, $\langle \mathcal{I}|\mathcal{G}\rangle {\rm and}\langle \mathcal{I}|\mathcal{I}\rangle$.

We time propagate the coefficients using an explicit fourth order Runge–Kutta method with adaptive step size. In order to absorb the wavefunction at the box boundaries we use infinite range exterior complex scaling (irECS) [23]. We typically choose simulation boxes larger than the spatial extent of the Hartree–Fock orbitals and start absorption after the Hartree–Fock orbitals vanish. This implies that it suffices to complex scale only one of the N coordinates.

The cost of time propagation scales with the number of ionic states (say n_I) as n_I² and it is independent of the exact number of electrons. This makes basis sets of the kind (5) attractive for treating many electron systems.

In order to solve the TDSE we need to evaluate various operators. Firstly, we introduce several generalized reduced density matrices with the help of creation ($a_{k}^{\dagger }$) and annihilation (a_k) operators on the single particle state ${{\phi }_{k}}$. A pth order generalized reduced density matrix between the (N-1) particle ionic functions is given by:

$$\begin{eqnarray}&&\rho _{{{k}_{1}},\ldots ,{{k}_{p}},{{l}_{1}},\ldots ,{{l}_{p}}}^{IJ}=\langle I|a_{{{k}_{1}}}^{\dagger }...a_{{{k}_{p}}}^{\dagger }{{a}_{{{l}_{1}}}}...{{a}_{{{l}_{p}}}}|J\rangle .\end{eqnarray} \tag{ 12 }$$

Similarly, we define generalized Dyson coefficients between the N-particle neutral wavefunctions and (N-1) particle ionic wavefunctions as

$$\begin{eqnarray}&&\eta _{{{k}_{1}},\ldots ,{{k}_{p}},{{l}_{1}},\ldots ,{{l}_{p-1}}}^{\mathcal{G}J}=\langle \mathcal{G}|a_{{{k}_{1}}}^{\dagger }...a_{{{k}_{p}}}^{\dagger }{{a}_{{{l}_{1}}}}...{{a}_{{{l}_{p-1}}}}|J\rangle .\end{eqnarray} \tag{ 13 }$$

With the help of these objects, we can present the final form of the matrix elements. The overlap matrix blocks have the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle \mathcal{G}|\mathcal{G}\rangle & = & 1 \\ \langle \mathcal{G}|\mathcal{I}\rangle & = & \eta _{k}^{\mathcal{G}I}\langle {{\phi }_{k}}|i\rangle \\ \langle \mathcal{I}|\mathcal{J}\rangle & = & \langle i|j\rangle \langle I|J\rangle -\langle i|{{\phi }_{l}}\rangle \rho _{kl}^{IJ}\langle {{\phi }_{k}}|j\rangle \\ {} & = & \langle i|j\rangle {{\delta }_{IJ}}-\langle i|{{\phi }_{l}}\rangle \rho _{kl}^{IJ}\langle {{\phi }_{k}}|j\rangle , \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where $\eta _{k}^{\mathcal{G}I}$ can be identified with the Dyson orbital coefficients with respect to the Hartree–Fock orbitals and $\rho _{kl}^{IJ}$ are the one particle reduced density matrices.

Any exchange-symmetric single particle operator can be written as

$$\begin{eqnarray}&&\hat{T}=\hat{t}(1)+\hat{t}(2)+...+\hat{t}(N),\end{eqnarray} \tag{ 15 }$$

where $\hat{t}(u)$ is the single particle operator corresponding to the coordinate u. Matrix elements of $\hat{T}$ are

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle \mathcal{G}|\hat{T}|\mathcal{G}\rangle & = & \rho _{kl}^{\mathcal{G}\mathcal{G}}\langle {{\phi }_{k}}|\hat{t}|{{\phi }_{l}}\rangle \\ \langle \mathcal{G}|\hat{T}|\mathcal{I}\rangle & = & \eta _{k}^{\mathcal{G}I}\langle {{\phi }_{k}}|\hat{t}|i\rangle +\eta _{klm}^{\mathcal{G}I}\langle {{\phi }_{k}}|\hat{t}|{{\phi }_{m}}\rangle \langle {{\phi }_{l}}|i\rangle \\ \langle \mathcal{I}|\hat{T}|\mathcal{J}\rangle & = & \langle i|j\rangle \langle I|\hat{t}|J\rangle +\langle i|\hat{t}|j\rangle \langle I|J\rangle -\langle i|\hat{t}|{{\phi }_{l}}\rangle \rho _{kl}^{IJ}\langle {{\phi }_{k}}|j\rangle \\ {} & {} & -\langle i|{{\phi }_{l}}\rangle \rho _{kl}^{IJ}\langle {{\phi }_{k}}|\hat{t}|j\rangle -\langle i|{{\phi }_{c}}\rangle \langle {{\phi }_{a}}|j\rangle \langle {{\phi }_{b}}|\hat{t}|{{\phi }_{d}}\rangle \rho _{abcd}^{IJ}, \\ \end{array}\end{eqnarray} \tag{ 16 }$$

where $\eta _{klm}^{\mathcal{G}I}$ are the three index generalized Dyson coefficients, equation (13), and $\rho _{abcd}^{IJ}$ are the two particle reduced density matrices, equation (12).

Finally, the two particle operators

$$\begin{eqnarray}&&{{\hat{V}}^{(2)}}=\mathop{\sum }\limits_{i\lt j}\hat{v}(ij)\end{eqnarray} \tag{ 17 }$$

have the matrix elements

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} \langle \mathcal{G}|{{{\hat{V}}}^{(2)}}|\mathcal{G}\rangle =\frac{1}{2}\rho _{abcd}^{\mathcal{G}\mathcal{G}}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{c}}{{\phi }_{d}}\rangle \\ \langle \mathcal{G}|{{{\hat{V}}}^{(2)}}|\mathcal{I}\rangle =\eta _{klm}^{\mathcal{G}I}\langle {{\phi }_{k}}{{\phi }_{l}}|\hat{v}|{{\phi }_{m}}i\rangle +\frac{1}{2}\eta _{abcde}^{\mathcal{G}I}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{d}}{{\phi }_{e}}\rangle \langle {{\phi }_{c}}|i\rangle \\ \langle \mathcal{I}|{{{\hat{V}}}^{(2)}}|\mathcal{J}\rangle =\frac{1}{2}\rho _{abcd}^{\mathcal{I}\mathcal{J}}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{c}}{{\phi }_{d}}\rangle \langle i|j\rangle +\mathop{\underbrace{\rho _{kl}^{IJ}\langle {{\phi }_{k}}i|\hat{v}|{{\phi }_{l}}j\rangle }}\limits_{{\rm Direct}\;{\rm term}}-\mathop{\underbrace{\rho _{kl}^{IJ}\langle {{\phi }_{k}}i|\hat{v}|j{{\phi }_{l}}\rangle }}\limits_{{\rm Standard}\;{\rm exchange}\;{\rm term}} \\ -\mathop{\underbrace{\rho _{abcd}^{IJ}\langle {{\phi }_{a}}i|\hat{v}|{{\phi }_{c}}{{\phi }_{d}}\rangle \langle {{\phi }_{b}}|j\rangle -\rho _{abcd}^{IJ}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{c}}j\rangle \langle i|{{\phi }_{d}}\rangle -\frac{1}{2}\rho _{abcdef}^{IJ}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{d}}{{\phi }_{e}}\rangle \langle i|{{\phi }_{f}}\rangle \langle {{\phi }_{c}}|j\rangle }}\limits_{{\rm Other}\;{\rm exchange}\;{\rm terms}\;{\rm due}\;{\rm to}\;{\rm non}-{\rm orthogonality}}, \\ \end{array}\end{eqnarray*} \tag{ 18 }$$

where $\eta _{abcde}^{\mathcal{G}I}$ are the five index generalized Dyson coefficients, equation(13) and $\rho _{abcdef}^{IJ}$ are the three particle reduced density matrices, equation (12). Although it appears from equation (18) that the necessity of $\rho _{abcdef}^{IJ}$, $\eta _{abcde}^{\mathcal{G}I}$ leads to large memory requirements, we must point out that the contractions $\rho _{abcdef}^{IJ}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{d}}{{\phi }_{e}}\rangle$ and $\eta _{abcde}^{\mathcal{G}I}\langle {{\phi }_{a}}{{\phi }_{b}}|\hat{v}|{{\phi }_{d}}{{\phi }_{e}}\rangle$ can be made while computing $\rho _{abcdef}^{IJ}$ and $\eta _{abcde}^{\mathcal{G}I}$ itself, thereby storing only simple matrices and vectors.

In order to compute the two-electron integrals, we first project the Hartree–Fock orbitals onto a single center expansion:

$$\begin{eqnarray}&&{{\phi }_{k}}(\vec{r})=\mathop{\sum }\limits_{{{q}_{k}}{{l}_{k}}{{m}_{k}}}{{d}_{{{q}_{k}}{{l}_{k}}{{m}_{k}}}}{{Y}_{{{l}_{k}}{{m}_{k}}}}(\theta ,\phi ),\end{eqnarray} \tag{ 19 }$$

where q_k are radial quadrature points, ${{l}_{k}},{{m}_{k}}$ refer to the angular momentum functions and use these expansions with the multi-pole expansion:

$$\begin{eqnarray}&&\frac{1}{|{{{\vec{r}}}_{1}}-{{{\vec{r}}}_{2}}|}=\mathop{\sum }\limits_{LM}\frac{4\pi }{2L+1}\frac{r_{\lt }^{L}}{r_{\gt }^{L+1}}{{Y}_{LM}}({{\theta }_{1}},{{\phi }_{1}})Y_{LM}^{*}({{\theta }_{2}},{{\phi }_{2}})\end{eqnarray} \tag{ 20 }$$

with ${{r}_{\lt }}={\rm min} ({{r}_{1}},{{r}_{2}})$ and ${{r}_{\gt }}={\rm max} ({{r}_{1}},{{r}_{2}})$.The limits for the multi-pole expansion are determined by the angular momenta in the one-electron numerical basis and the single center expansion for the molecular orbitals ${{\phi }_{k}}$. No other truncation schemes are employed. These two particle operators pose a challenge for efficient computation. While the direct term is relatively easy to handle, the exchange terms consume a major portion of the Hamiltonian setup time.

The overlap matrix (14) is not a standard finite element overlap matrix, which would be banded and would allow for efficient application of the inverse. Rather, non-orthogonality between the active electron basis and the Hartree–Fock orbitals leads to extra cross terms that destroy the banded structure in general and complicate the computation of the inverse. However, the inverse of the overlap can still be computed efficiently using low rank updates. We use here the Woodbury formula [25], according to which the inverse of a modified matrix of the form $({{S}_{0}}-U\Lambda {{U}^{\dagger }})$ can be computed as:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{S}^{-1}} & = & {{({{S}_{0}}-U\Lambda {{U}^{\dagger }})}^{-1}} \\ {} & = & S_{0}^{-1}-S_{0}^{-1}U{{({{U}^{\dagger }}S_{0}^{-1}U-{{\Lambda }^{-1}})}^{-1}}{{U}^{\dagger }}S_{0}^{-1}. \\ \end{array}\end{eqnarray} \tag{ 21 }$$

As an example, with 2 ionic states and 1 neutral the overlap matrix (14) can be cast in the form:

$$\begin{eqnarray}S=\mathop{\underbrace{\left( \begin{array}{ccccccccccccccc} {{s}_{0}} & 0 & 0 \\ 0 & {{s}_{0}} & 0 \\ 0 & 0 & 1 \\ \end{array} \right)}}\limits_{{{S}_{0}}}-\mathop{\underbrace{\left( \begin{array}{ccccccccccccccc} u & 0 & 0 \\ 0 & u & 0 \\ 0 & 0 & 1 \\ \end{array} \right)}}\limits_{U}\mathop{\underbrace{\left( \begin{array}{ccccccccccccccc} {{\rho }^{11}} & {{\rho }^{12}} & {{\eta }^{\mathcal{G}1}} \\ {{\rho }^{21}} & {{\rho }^{22}} & {{\eta }^{\mathcal{G}2}} \\ {{[{{\eta }^{\mathcal{G}1}}]}^{T}} & {{[{{\eta }^{\mathcal{G}2}}]}^{T}} & 0 \\ \end{array} \right)}}\limits_{\Lambda }{{U}^{\dagger }}\end{eqnarray} \tag{ 22 }$$

which is suitable for the Woodbury formula (21). Here, ${{({{s}_{0}})}_{ij}}=\langle i|j\rangle$ and ${{u}_{ik}}=\langle i|{{\phi }_{k}}\rangle$. Let n_a denote the number of active electron basis functions $|i\rangle$ and n_hf be the number of Hartree–Fock orbitals ${{\phi }_{k}}$ that is much smaller that n_a. The overlap matrix s₀ has dimensions ${{n}_{a}}\times {{n}_{a}}$ but is narrowly banded, and the dimensions of matrix u are ${{n}_{a}}\times {{n}_{hf}}$.

Let n_I be the number of ionic states. Then the overlap S and S₀ are $({{n}_{I}}{{n}_{a}}+1)\times ({{n}_{I}}{{n}_{a}}+1)$ matrices, where the inverse of S₀ can be easily applied. The matrix U is $({{n}_{I}}{{n}_{a}}+1)\times ({{n}_{I}}{{n}_{hf}}+1)$ and Λ is $({{n}_{I}}{{n}_{hf}}+1)\times ({{n}_{I}}{{n}_{hf}}+1)$. This low rank structure of the correction terms can be utilized to compute the inverse efficiently by using the Woodbury formula.

Anti-symmetrization and non-orthogonality of the active electron basis with respect to the Hartree–Fock orbitals may render our basis linearly dependent. If the $\{|i\rangle \}$-basis is near-complete w.r.t. the HF-orbital basis

$$\begin{eqnarray}&&\langle {{\phi }_{k}}|i\rangle \ {{[s_{0}^{-1}]}_{ij}}\ \langle j|{{\phi }_{l}}\rangle \approx {{\delta }_{kl}},\end{eqnarray} \tag{ 23 }$$

it is possible to find coefficients ${{c}_{i,I}}$ such that

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i,I}{{c}_{i,I}}\mathcal{A}[|i\rangle |I\rangle ]\approx 0.\end{eqnarray} \tag{ 24 }$$

An obvious case where this happens is when one and the same HF orbital ${{\phi }_{{{k}_{0}}}}$ appears in all the ionic determinants. For a linear combination ${{\sum }_{i}}{{c}_{i,I}}|i\rangle \approx {{\phi }_{{{k}_{0}}}}$ anti-symmetrization renders equation (24) near zero. As a result, the overlap matrix becomes non-invertible. A possible solution would be to orthogonalize the active electron basis with respect to the Hartree–Fock orbitals. But this is not an easily implementable solution with a CI ionic basis. For each determinant, the set of Hartree–Fock orbitals with respect to which the active electron basis must be orthogonal is different.

As an alternative solution we use a generalization of the Woodbury formula (21) to compute the inverse of a matrix only on the subspace of the non-zero eigenvectors of the matrix. Let Z denote the ${{n}_{0}}\times {{n}_{z}}$ matrix of eigenvectors with near-zero eigenvalues ${{z}_{p}}\lt \epsilon ,p=1,\ldots ,{{n}_{z}}$ of the generalized eigenvalue problem

$$\begin{eqnarray}&&SZ={{S}_{0}}Z{{d}_{z}},\end{eqnarray} \tag{ 25 }$$

with d_z denoting the diagonal matrix of the eigenvalues z_p and Z satisfying the orthonormality relation ${{Z}^{\dagger }}{{S}_{0}}Z={\bf 1}$. In general, there will be comparatively few such eigenvectors ${{n}_{z}}\ll {{n}_{0}}$ and these can be easily determined by an iterative solver. We can remove these singular vectors from our calculation by the projector

$$\begin{eqnarray}&&Q=1-Z{{Z}^{\dagger }}{{S}_{0}}.\end{eqnarray} \tag{ 26 }$$

The projector property ${{Q}^{2}}=Q$ can be easily verified. As the projector refers to the generalized eigenvalue problem with ${{S}_{0}}\ne {\bf 1}$, Q is not an orthogonal projector, that is ${{Q}^{\dagger }}\ne Q$. We define a pseudo-inverse ${{\tilde{S}}^{-1}}$ of S on the subspace of generalized eigenvectors with non-zero eigenvalues with the property

$$\begin{eqnarray}&&{{\tilde{S}}^{-1}}SQ=Q.\end{eqnarray} \tag{ 27 }$$

One can verify directly that the generalized Woodbury formula

$$\begin{eqnarray}&&{{\tilde{S}}^{-1}}=QS_{0}^{-1}\left[ 1-U{{({{U}^{\dagger }}QS_{0}^{-1}U-\Lambda )}^{-1}}{{U}^{\dagger }}QS_{0}^{-1} \right]\end{eqnarray} \tag{ 28 }$$

satisfies the definition (27). The matrix $({{U}^{\dagger }}QS_{0}^{-1}U-\Lambda )$ is invertible on all vectors appearing in ${{U}^{\dagger }}Q$ to its right, as exactly the singular vectors are removed by the projector Q. Apart from the necessity to determine Z during setup, the correction does not significantly increase the operations count for the inverse overlap.

In [22], we had shown that when an electron is treated with a restricted basis, for example, in terms of a few bound states, the length gauge is a more natural gauge. Compared to pure velocity gauge, the coupled channel computations converge quickly in mixed gauge with length gauge spanning the region of the ionic states and velocity gauge thereafter for asymptotics. In this current work, we use continuous gauge switching, detailed in [22], for its easy implementation. Starting from the length-gauge, we solve the TDSE after applying the gauge transformation

$$\begin{eqnarray}{{U}_{c}}=\left\{ \begin{array}{ccccccccccccccc} 1 & {\rm for}\;r\leqslant {{r}_{{\rm g}}}, \\ {\rm exp} \left[ {\rm i}\vec{A}(t)\cdot \mathop{\sum }\limits_{j=1}^{N}\widehat{{{r}_{j}}}\left( {{r}_{j}}-{{r}_{{\rm g}}} \right) \right] & {\rm for}\;r\gt {{r}_{{\rm g}}}. \\ \end{array} \right.\end{eqnarray} \tag{ 29 }$$

Here, r_g is the gauge radius that separates the length gauge and velocity gauge regions.

The computation of photoelectron spectra is expensive for two reasons. (1) The asymptotic part of the wavefunction is needed to extract photoelectron spectra, which means large simulation boxes to preserve the asymptotic part and to avoid any numerical reflections that may corrupt the wavefunction. (2) Single continuum states are needed into which the wavefunction must be decomposed, in order to obtain photoelectron spectra. These two problems are circumvented in a recently developed method tSURFF [19, 20]: one computes the spectra by a time integration of electron flux through a surface defined by a radius R_c called the tSURFF radius. The Coulomb potential is smoothly turned off before R_c, which implies that the scattering solutions thereafter are the well known Volkov solutions. R_c becomes a convergence parameter, and by varying this radius, one can compute spectra to a given accuracy. This method has been explained in detail in previous works for single ionization in [19] and for double ionization in [20]. A proposal for extension of this method for single ionization of multi-electron systems has been outlined in [20]. We describe here the application of the method with our coupled channels setup.

Let ${{\chi }_{k}}$ be the scattering solutions which take the form of Volkov solutions beyond R_c and $\Psi (T)$ be the wavefunction at some large time T. According to tSURFF for single electron systems, photoelectron spectra can be computed as ${{\sigma }_{k}}=|{{b}_{k}}{{|}^{2}}$ with b_k defined as:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{b}_{k}} & = & \langle {{\chi }_{k}}(T)|\Theta ({{R}_{c}})|\Psi (T)\rangle \\ {} & = & {\rm i}\int _{0}^{T}{\rm d}t\langle {{\chi }_{k}}(t)|\left[ -\frac{1}{2}\vartriangle +{\rm i}\vec{A}(t)\cdot \vec{\nabla },\Theta ({{R}_{c}}) \right]|\Psi (t)\rangle . \\ \end{array}\end{eqnarray} \tag{ 30 }$$

Here $\Theta ({{R}_{c}})$ is a Heaviside function that is unity for $r\gt {{R}_{c}}$ and 0 elsewhere.

This formulation can be easily extended to the N electron problem in a coupled channels setup. In the present setup, we mostly take a set of ionic bound states for the ionic basis. These states have a finite extent. We may choose R_c such that the electrons described by the ionic basis vanish by R_c which means all the exchange terms in the Hamiltonian vanish after R_c. The remaining direct potential can be smoothly turned off just as the Coulomb potential. This implies that the wavefunction beyond R_c evolves by the Hamiltonian:

$$\begin{eqnarray}&&H(r\gt {{R}_{c}})={{H}_{{\rm ion}}}\otimes \hat{1}+{{\hat{1}}_{{\rm ion}}}\otimes \left[ -\frac{1}{2}\vartriangle +{\rm i}\vec{A}(t)\cdot \vec{\nabla } \right]\end{eqnarray} \tag{ 31 }$$

that allows for a complete set of solutions of the form:

$$\begin{eqnarray}&&{{\xi }_{c,k}}({{\vec{r}}_{1}},\ldots ,{{\vec{r}}_{n}},t)=\mathcal{A}\left[ {{\kappa }_{c}}({{{\vec{r}}}_{1}},\ldots ,{{{\vec{r}}}_{n-1}},t)\otimes {{\chi }_{k}}({{{\vec{r}}}_{n}},t) \right],\end{eqnarray} \tag{ 32 }$$

where H_ion is the single ionic Hamiltonian and ${{\kappa }_{c}}(t)$ are time dependent ionic channel functions solving the TDSE

$$\begin{eqnarray}&&{\rm i}\frac{\partial {{\kappa }_{c}}(t)}{\partial t}={{\hat{H}}_{{\rm ion}}}{{\kappa }_{c}}(t)\end{eqnarray} \tag{ 33 }$$

within the ansatz in terms of field-free ionic states

$$\begin{eqnarray}&&|{{\kappa }_{c}}(t)\rangle =\mathop{\sum }\limits_{I}|I\rangle {{d}_{cI}}(t).\end{eqnarray} \tag{ 34 }$$

With the help of the ${{\xi }_{c,k}}$, channel resolved photoelectron spectra can be computed as

$$\begin{eqnarray}&&{{\sigma }_{c,k}}=|\langle {{\xi }_{c,k}}({{\vec{r}}_{1}},\ldots ,{{\vec{r}}_{n}},T)|\Theta ({{R}_{c}})|\Psi ({{\vec{r}}_{1}},\ldots ,{{\vec{r}}_{n}},T)\rangle {{|}^{2}}\end{eqnarray} \tag{ 35 }$$

and the asymptotic decomposition of Ψ in terms of ${{\xi }_{c,k}}$ is obtained as

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \langle {{\xi }_{c,k}}({{{\vec{r}}}_{1}},\ldots ,{{{\vec{r}}}_{n}},T)|\Theta ({{R}_{c}})|\Psi ({{{\vec{r}}}_{1}},\ldots ,{{{\vec{r}}}_{n}},T)\rangle \\ \qquad ={\rm i}\int _{0}^{T}{\rm d}t\langle \mathcal{A}\left[ {{\kappa }_{c}}({{{\vec{r}}}_{1}},\ldots ,{{{\vec{r}}}_{n-1}},t)\otimes {{\chi }_{k}}({{{\vec{r}}}_{n}},t) \right]|\left[ -\frac{1}{2}{{\vartriangle }_{n}}+{\rm i}\vec{A}(t)\cdot {{{\vec{\nabla }}}_{n}},{{\Theta }_{n}}({{R}_{c}}) \right]|\Psi ({{{\vec{r}}}_{1}},\ldots ,{{{\vec{r}}}_{n}},t)\rangle \\ \qquad ={\rm i}\int _{0}^{T}{\rm d}t\langle {{\chi }_{k}}({{{\vec{r}}}_{n}},t)|\left[ -\frac{1}{2}{{\vartriangle }_{n}}+{\rm i}\vec{A}(t)\cdot {{{\vec{\nabla }}}_{n}},{{\Theta }_{n}}({{R}_{c}}) \right]|{{\zeta }_{c}}({{{\vec{r}}}_{n}},t)\rangle , \\ \end{array}\end{eqnarray} \tag{ 36 }$$

where we introduced the time-dependent Dyson orbitals

$$\begin{eqnarray}&&{{\zeta }_{c}}({{\vec{r}}_{n}},t)\;:=\;\langle {{\kappa }_{c}}({{\vec{r}}_{1}},\ldots ,{{\vec{r}}_{n-1}},t)|\Psi ({{\vec{r}}_{1}},\ldots ,{{\vec{r}}_{n}},t)\rangle .\end{eqnarray} \tag{ 37 }$$

The commutator of the derivatives with the Heaviside function Θ gives δ-like terms involving values and derivatives of Ψ at the surface $|\vec{r}|={{R}_{c}}$. As we choose R_c such that the Hartree–Fock orbitals vanish by then, we do not need to consider the exchange terms in computing ${{\zeta }_{c}}$. Along with time propagating the N electron problem, one needs to also time propagate the ionic problem (33).

A detailed discussion of performance and intrinsic limitations of the tSURFF method is contained in [19, 20]. We here summarize the main points of this discussion. The strength of tSURFF lies (a) in a dramatic reduction of the required numerical box sizes to compute accurate spectra and (b) in the fact that no scattering states are needed for spectral analysis. As the asymptotic scattering information is generated during time-propagation rather than by solving an independent stationary problem, propagation times must be long enough for all relevant processes to terminate and for all electrons to pass through the surface. This favors the application of the method for fast processes. For slow processes like emission at near zero electron energy or the decay of long-lived resonances purely stationary methods or methods that combine solutions of the TDSE during pulse with a stationary analysis after the end of the pulse may become advantageous. Also, if small boxes are used, the capability for representing very extended objects like Rydberg states is limited by the box size. We will illustrate these points below when discussing photo-emission from the helium atom at short wave length.

As we solve the non-relativistic TDSE, the total spin of the system is conserved during the time evolution. We can therefore remove the spin degree of freedom through suitable linear combinations of the anti-symmetrized products in the basis (5) to enforce a particular spin symmetry. This reduces the size of our basis. We consider only singlet spin symmetric systems in this work. As an example, we show how singlet spin symmetry can be enforced. Let $\uparrow$ and $\downarrow$ indicate the spin states $\pm \frac{1}{2}$ associated with a spatial function. Choosing linear combinations of the kind:

$$\begin{eqnarray}&&\mathcal{A}[|I\rangle |i\rangle ]:=\frac{\mathcal{A}\left[ |{{I}^{\uparrow }}\rangle |{{i}^{\downarrow }}\rangle \right]-\mathcal{A}\left[ |{{I}^{\downarrow }}\rangle |{{i}^{\uparrow }}\rangle \right]}{\sqrt{2}}\end{eqnarray} \tag{ 38 }$$

enforces singlet symmetry. This can be extended to creating linear combinations that enforce an arbitrary spin symmetry.

Helium is the largest atom that can be numerically treated in full dimensionality. With linearly polarized laser fields, the symmetry of the system can be used to reduce the problem to five-dimensions. The energies of the helium ionic states are $-2/{{n}^{2}}$ for principal quantum number n. The first two ionic states are separated by 1.5 a.u. in energy, which is large, for example, compared to a photon energy of 0.456 a.u. at 100 nm. This has been a motivation to treat helium as an effective single electron system with XUV and longer wavelengths in some earlier works, for example in [26]. We examine below the validity of treating helium as an effective single electron system, by comparing haCC calculations with full dimensional calculations at different wavelengths.

Figure 1 shows photoelectron spectra from helium with a 21 nm (ω = 2.174 a.u.), three cycle laser pulse with a peak intensity of ${{10}^{15}}\;{\rm W}\;{\rm c}{{{\rm m}}^{-2}}$. The one and two photon ionization peaks of 1s and 2p_z channel spectra are shown. The relative errors of haCC calculations are computed with respect to the full dimensional calculation. The single photon peak of the 1s channel is computed to a few percent accuracy, except for a feature around 1.3 a.u., with a single ionic state. The resonant feature can be identified with the 2s2p doubly excited state [27], which is reproduced to few percent accuracy with the addition of 2nd shell ionic states. While the position of the resonance is reproduced accurately in the calculations presented here, the propagation time was well below the life-time of this resonance which is reflected in the width of the feature that is well above the natural line width.

Zoom In Zoom Out Reset image size

The two photon peak of the 1s channel and the 2p_z channel spectra (figure 1) need more than a single ionic state and they could be computed only up to 15% accuracy even after inclusion of 9 ionic states ($n\leqslant 3$).

A broadband (few cycle) XUV pulse tends to excite the initial state into a band of final states which may include many correlated intermediate states. Here, the intrinsic limitations of any coupled channels approach that is based on ionic bound states only are exposed. Firstly, a correlated intermediate state with a bound character needs large number of ionic states to be correctly represented. Secondly, the ionic bound states based on Gaussian basis sets do not have the exact asymptotic behavior. This can lead to an inaccuracy in length gauge dipole matrix elements. Finally, the absence of ionic continuum states in our approach is another possible source of inaccuracy. Due to these limitations, we do not expect the shake-up channel spectra to be more accurate than 10–15%.

For obtaining long-lived resonance structures by a time-dependent method one must, as a general feature of such methods, propagate for times at least on the scale of the life time of the resonance. The only alternative is to independently solve the stationary resonant scattering problem and decompose the time-dependent solution after the end of the pulse into the corresponding scattering continuum. Solving the scattering problem, however, is a computationally very demanding task by itself. For obtaining the resonances with tSURFF, one can simply propagate until the resonance has decayed completely and all flux has passed through the surface where the flux is collected. At this point it should be remarked that the relevant information about resonances may be generated more efficiently by stationary methods like time-independent complex scaling [27, 28]. As the comparison in figure 1 is with the two-electron code where long propagation times become rather costly, we compare the spectra at time T = 60 laser cycles, where the resonances have not emerged yet.

With haCC, due to its very compact representation, we can easily propagated much longer to obtain the resonances to any desired accuracy. Figure 2 shows how the $n\leqslant 3$ resonances emerge with increasing propagation times in the 1s and 2p_z channels. For example, at the 2s2p decay width of $\Gamma =1.37\times {{10}^{-3}}$, one expects $76.5\%$ of the Auger electrons to have passed through the surface $|\vec{r}|=45$ at time T = 400, which increases to $95.2\%$ at T = 800. The ratio of 1.24 between these numbers closely matches the increase of mass in the 2s2p peak by 1.23 in figure 2. The positions of the resonances are accurate on the level of 10⁻³ a.u. in energy [27], showing that the correlated states are well represented by the method. Here we had chosen a box size of 45 a.u. which is sufficient to represent doubly excited states from 2s2p to 2s6p states. Hence, only these states are seen in the spectra. One can obtain the higher excited states by increasing the box size, at the penalty of a larger discretization and somewhat longer propagation times, see discussion in section 2.7. Similar as long lived resonances, threshold behavior of the spectrum near energy zero only emerges with propagation time, figure 2. Further distortions near threshold are due to the effective truncation of the Coulomb tail in the absorbing region. In the present example, these effects produce an artefact at energies $\lesssim 0.02$ a.u. in the right panel of figure 2. Accuracies at the lowest energies can be pushed by increasing both, simulation box size and propagation times.

Zoom In Zoom Out Reset image size

Figure 3 shows total photoelectron spectra from helium at 200 and 400 nm wavelengths. The exact laser parameters are indicated in the figure captions. At 200 nm (ω = 0.228 a.u.), the ionization threshold is four photons. A single ionic state calculation produces spectra that are 10% accurate with respect to a full dimensional calculation. Addition of second and third shell ionic states improves the accuracy of the spectra to few percent level in the important regions of the spectrum. At 400 nm (ω = 0.114 a.u.), the ionization threshold is eight photons. Also here, a single ionic state computation produces spectra that are 10% accurate with respect to a full dimensional calculation. Addition of more ionic states, does not improve the accuracy further. This is possibly due to the missing continuum of the second electron that is needed to fully describe the polarization of the ionic core.

Zoom In Zoom Out Reset image size

At longer wavelengths, we find that single ionic state computations are sufficient to produce spectra accurate on the level of 10%. This is consistent with the knowledge that at longer wavelengths, it is the ionization thresholds that play a more important role in determining the ionization yields compared to the exact electronic structure. Our findings show that helium at long wavelengths can be approximated as a single channel system.

Beryllium is a four electron system that is often treated as a two electron system due to the strong binding of its inner two electrons. The third ionization potential of beryllium is 153.8961 eV [29]. With photon energies below this third ionization potential, it can be safely treated as an effective two electron system. This allows us to have a benchmark for our spectra by adapting the simple Coulomb potential to an effective potential in our two electron code. We use the effective potential given in [6] for our benchmark calculations. We refer to these as 'effective-2e' calculations.

Table 1 lists the energies of the first 8 ionic states of beryllium relative to the ground ionic state. As the ionic excitation energies are much smaller than in Helium one would expect inter-channel couplings to play a greater role.

Table 1.  Energies of the used single ionic states of beryllium relative to the ground state ion. The COLUMBUS [21] ionic states are computed at MR-CISD level with aug-cc-pvtz basis.

Ionic state	NIST database (eV)	Columbus energies (eV)
$1{{{\rm s}}^{2}}2{\rm s}$	0.0	0.0
$1{{{\rm s}}^{2}}2{\rm p}$	3.9586	3.9767
$1{{{\rm s}}^{2}}3{\rm s}$	10.9393	10.9851
$1{{{\rm s}}^{2}}3{\rm p}$	11.9638	12.1407

Figure 4 shows photoelectron spectra from beryllium with 21 and 200 nm wavelength laser pulses. The exact parameters are indicated in the figure caption. The relative errors of spectra from the haCC calculations are computed with respect to the effective-2e calculations.

Zoom In Zoom Out Reset image size

At 21 nm, the one and two photon ionization peaks of ground state channel spectra are shown. Here, the single photon ionization process itself needs more than the ground ionic state to produce accurate photoelectron spectra. Adding more ionic states improves the accuracy to a few percent level. We find that the close energetic spacing of beryllium ionic states leads to a greater possibility of inter-channel coupling, which requires more than the ground ionic state to be well represented.

Also, at 200 nm we need more than the ground ionic state to compute realistic spectra. With $1{{{\rm s}}^{2}}2{\rm s}$, $1{{{\rm s}}^{2}}2{\rm p}$ ionic states included, the spectra produced have 20% accuracy with respect to the benchmark calculation. With the addition of $1{{{\rm s}}^{2}}3{\rm s}$ and $1{{{\rm s}}^{2}}3{\rm p}$ states, a structure similar to the one predicted by the benchmark calculation develops around 10 eV. This structure may be identified with the lowest resonance $1{{{\rm s}}^{2}}2{\rm p}3{\rm s}$ at 10.71 eV [6]. The coupled channels calculations with the number of ionic states considered here, however do not reproduce the structure on the second peak exactly. This points to a feature of a coupled channels basis that the correct representation of a strongly correlated state that has bound character requires a large number of ionic states. As an alternative strategy, one can explicitly include the correlated state of importance into the basis, if it can be pre-computed, on the same footing as the correlated ground state.

It has been shown through examples in section 4.1 that helium can be modeled as single channel system at longer wavelengths. Lithium, the smallest alkali metal, also has been successfully modeled as a single electron system in an effective potential, for example in [30]. We find that beryllium needs at least two ionic states, $1{{{\rm s}}^{2}}2{\rm s}$ and $1{{{\rm s}}^{2}}2{\rm p}$, for a realistic modeling. It serves as a first simple example where single electron models break down and multiple channels need to be considered.

In figure 5, we present the single photon ionization cross-sections as a function of photon energy from our haCC method and compare them with the cross-sections calculated with TD-RASCI method [11] and R-matrix method [31] and with experimental results from [32]. The cross-sections in our time dependent approach are computed using the equation (51) given in [18] with which the N photon ionization cross-section, ${{\sigma }^{(N)}}$ in units ${\rm c}{{{\rm m}}^{2N}}/{{{\rm s}}^{N-1}}$ can be computed as:

$$\begin{eqnarray}&&{{\sigma }^{(N)}}={{(8\pi \alpha )}^{N}}{{\left( \frac{3.5\times {{10}^{16}}}{I} \right)}^{N}}{{\omega }^{N}}\Gamma a_{0}^{2N}t_{0}^{N-1},\end{eqnarray} \tag{ 42 }$$

where I is the intensity in ${\rm W}\;{\rm c}{{{\rm m}}^{-2}}$, ω is the laser frequency in a.u., α is the fine structure constant and a₀, t₀ are atomic units of length and time respectively. Γ is the total ionization rate which is computed in a time dependent approach by monitoring the rate at which the norm of the wavefunction in a certain inner region drops. As we are computing the rate, the exact size of the inner region does not play a role. The norm drop reaches a steady state irrespective of the inner region size. We used for our computations presented here a 40-cycle continuous wave laser pulse with a three cycle ${{{\rm cos} }^{2}}$ ramp up and ramp down and with an intensity of ${{10}^{12}}\;{\rm W}\;{\rm c}{{{\rm m}}^{-2}}$. We checked convergence with respect to the pulse duration and the inner region size, and the computations are converged well below the differences seen by addition of ionic states in the basis.

Zoom In Zoom Out Reset image size

All the theoretical methods agree with each other qualitatively, though there are differences on the level of 5–10% quantitatively. The experimental results from [32] have error bars on the level ≲10% (0.1 Mb) which are not shown here. All the theoretical results lie in this range except at low energies. In the higher photon energy range, 30–60 eV, the haCC results and the R-matrix results are in good agreement compared to the TD-RASCI. In the haCC calculations, including more than 4 ionic states does not change the cross-sections. In the photon energy range 20–30 eV, the haCC computations with 5 and 8 ionic states are in good agreement with TD-RASCI results compared to the R-matrix results. In this energy range, the cross-sections from haCC calculations show a dependence on the number of ionic states included. This modulation may be attributed to the presence of auto-ionizing states in this region. Table 3 in [6] presents a list of resonances that appear in beryllium electronic structure. The first ionization potential is 9.3 eV. With photon energies around 20 eV, the resulting photoelectron reaches continuum region where a number of resonances are present. As correlated resonances need many ionic states to be well represented in a coupled channels basis, this may explain the dependence of the cross-section on the number of ionic states in 20–30 eV photon range.

The hydrogen molecule in linearly polarized laser fields parallel to the molecular axis, with fixed nuclei has the same symmetry as helium in linearly polarized laser fields. The off-centered nuclear potential, however, increases the angular momenta requirement when treated with a single center expansion. While the number of basis functions can be reduced through a choice of a more natural coordinate system like prolate spheroidal coordinates for diatomics [33], the challenge of computing two electron integrals remains. In the case of hydrogen molecule at equilibrium internuclear distance (R₀ = 1.4 a.u.), a calculation with single center expansion easily converges, as the proton charges do not significantly distort the spherical symmetry of the electron cloud. As a benchmark for spectra, we use results from a full dimensional calculation, that expands the wavefunction in a single center basis.

Figure 6 shows photoelectron spectra from H₂ at 21 nm wavelength. The exact laser parameters are given in the figure caption. The ground state ($1{{\sigma }_{{\rm g}}}$) and first excited state ($1{{\sigma }_{u}}$) channel spectra are shown. We find that, at this wavelength, a single ionic state is not sufficient to produce accurate photoelectron spectra. With the addition of more ionic states, there is a systematic improvement in the accuracy of the calculations. With 11 lowest σ and π ionic states included, we obtain an accuracy of about 10% for the $1{{\sigma }_{{\rm g}}}$ channel. The single photon ionization to the shake-up channel $1{{\sigma }_{u}}$ is also computed to a few percent accuracy with 11 ionic states. We find that the single ionization continuum of H₂ is more complex than helium and it needs more than a single ionic state.

Zoom In Zoom Out Reset image size

With 4 and 6 ionic states, we find artefacts on the two photon peaks. This is a result of a part of the COLUMBUS neutral ground state $|\mathcal{G}\rangle$ appearing in the eigenvalue spectrum of the Hamiltonian as a spurious doubly excited state $|s\rangle$. Let ${{\Pi }_{C}}$ be the projector onto the subspace spanned by the coupled-channels basis $\mathcal{A}[|i\rangle |I\rangle$. Then parts of the correlation contained in $|\mathcal{G}\rangle$ cannot be presented in that basis such that a non-zero correlated state

$$\begin{eqnarray}&&|s\rangle =(1-{{\Pi }_{C}})|\mathcal{G}\rangle \ne 0\end{eqnarray} \tag{ 43 }$$

appears at elevated energies. This spurious correlated state moves to higher energy with addition of ionic states. A straight forward solution to this problem is to compute this state and project it out from the basis. But this would require locating the spurious state in the eigenvalue spectrum, which is very demanding for large Hamiltonians. Fortunately, by their dependence on the number of ionic states, artefacts of this kind are easily detected and can be moved out of the region of interest by using sufficiently many ionic states. Such artefacts are a natural consequence of any ansatz of the kind (5) and need to be monitored.

Figure 7 shows total photoelectron spectra at 200 and 400 nm wavelengths. At 200 nm, spectra are accurate up to 10% with 2 ionic states. Addition of more ionic states helps reproduce additional resonant features in the spectrum. Also at 400 nm, 2 ionic states are sufficient to compute spectra that are accurate on 10% level, except for the resonant features. Inclusion of up to 6 ionic states reproduces the feature around 0.62 a.u. in the 400 nm spectrum, which may be attributed to second or third $^{1}\Sigma _{u}^{+}$ doubly excited state [34]. We find that with H₂ at longer wavelengths, ground ionic state is sufficient to compute realistic spectra and only for resonant features a large number ionic states is required.

Zoom In Zoom Out Reset image size

Figure 8 shows total ionization yield as a function of photon energy in the range 0.17–0.5 a.u. Results from haCC are compared with data available from other theoretical methods—time-dependent CI method from [35] and FNA-TDSE (fixed nuclei approximation) method from [7]. In addition, several points from our tSURFF-based full-2e method are included. The haCC calculations shown were performed using two ionic states, convergence was verified by performing four-state calculations at selected points. The vertical lines in the figure separate different multi-photon ionization regimes. The haCC, CI and full-2e are in fair agreement, while FNA-TDSE reproduces the threshold behavior, but severely, by up to an order of magnitude, deviates from the other calculations. The most conspicuous discrepancies between haCC and CI appear in the range $0.3\sim 0.4\;$ a.u. where CI exceeds haCC by about 20%. The discrepancies may be a result of the intrinsic limitations or the convergence of either calculation. For example, there are minor discrepancies in the ionization potential: the accurate value at H₂ equilibrium inter-nuclear distance (${{R}_{0}}=1.4\;$a.u.) is 0.6045 a.u.(table 1 in [36]),the ionization potential in [35] is 0.590 36 a.u., whereas in our calculations we obtain 0.6034 a.u. Also note that the results in [35] were shifted by 0.0092 a.u. in energy to match the resonance at 0.46 a.u. Although these differences are miniscule for energies they may indicate for somewhat larger deviations in the wave functions and the values of ionization potentials give a measure for the accuracy of the computations. Our full-2e computations that, in principle, could help to resolve the discrepancy are expensive and have not been pushed to an accuracy level which would allow to decide between the two results. However, we believe that the present level of agreement between haCC and CI is quite satisfactory and supports the validity of both approaches.

Zoom In Zoom Out Reset image size