Time-dependent convergent close coupling method for
molecular ionization in laser fields

Apparently, the only approach that allows one to correctly describe the nuclear and electronic dynamics during single ionization of a molecule without enormous expenditure of computational resources is the multi-configuration method.

To date, many variants of this approach have been proposed and used. The most advanced method appears to be the one developed by the group of F. Martin et al [1, 2]. This method is based on the use of pseudostates obtained by diagonalizing the Hamiltonian on a hybrid basis, assembled from functions of three types: Gaussian multicenter molecular orbitals, single-center Gaussian functions, and b-splines [2]. Multicenter Gaussian functions are used to approximate the wave function near the nuclear core. The advantages of using them are fast calculation of multicenter matrix elements and compatibility with programs for quantum chemical calculations that use the same basis. On the contrary, b-splines depending on the radial variable are used to approximate the far part of the wave function describing escaping electrons [3].

Nevertheless, despite all the advantages, this method also has disadvantages. The size of matrices and, accordingly, the speed of calculations for time evolution grows quadratically with the number of pseudostates. So increasing the energy resolution, or the maximum energy of emitted electrons, causes calculations to slow down dramatically. And for problems like Reconstruction of Attosecond Beating By Interference of Two-photon Transitions (RABBIT) with the resolution of nuclear states [4] or the observation of autoionization states, a high energy resolution is required.

Therefore, we developed a method that does not use pseudostates for the continuum.

Consider a two-electron molecule with the total electron spin $S$ and fixed nuclei. Time-dependent Schrd̈inger equation for the coordinate part of the wave function is

where the Hamiltonian can be written as

Here $\hat{h}_{0}(\mathbf{r})$ is the Hamiltonian of a singly ionized ion, $v(\mathbf{r}_{1},\mathbf{r}_{2})$ is the potential of electron-electron interaction, $\hat{w}(\mathbf{r},t)$ is the Hamiltonian of interaction with an external field.

Since we are interested in the ionization process, we will construct a multi-configuration wave function based on the function $\phi_{n}(\mathbf{r})$, that satisfy the stationary Schrd̈inger equation for the ion

The multi-configuration wave function can be written as

where $N_{s}$ is the number of ion states used. This structure of the function guarantees that after the emission of one electron, the remaining electron will end up in one of the ion’s own states.

If we project both sides of Eq. (1) onto $\phi_{n}(\mathbf{r}_{2})$, we obtain the equation

Substituting function (4) here, we obtain the equation

This equation is inconvenient because it contains two terms with time derivatives. To solve this problem, let’s introduce a new function

The introduction of a new function allows us to combine two terms with time derivatives into one.

To simplify the expression on the right side, consider the properties of the new function. The projections of the new function onto the basis functions satisfy the relation

from which it follows that it must satisfy the symmetry condition

A similar condition was once proposed to impose the convergent close-coupling [5], so in what follows we will call such relations the CCC symmetry.

Similar to the time-independent CCC [5], the same function $\Psi$ can be written in terms of different $\psi_{n}(\mathbf{r}_{1},t)$, so without loss of generality we can impose the CCC symmetry condition on the original function, i.e. put that $\langle\phi_{n}(\mathbf{r})|\psi_{k}(\mathbf{r},t)\rangle=(-1)^{S}\langle\phi_{k}(\mathbf{r})|\psi_{n}(\mathbf{r},t)\rangle$. Then $\langle\phi_{k}(\mathbf{r})|\bar{\psi}_{n}(\mathbf{r},t)\rangle=2\langle\phi_{k}(\mathbf{r})|\psi_{n}(\mathbf{r},t)\rangle$, and we can express the old function via the new one

Let us introduce another auxiliary function obtained by orthogonalization of the function with respect to the basis functions

With the use of the newly introduced functions and the CCC symmetry, we can rewrite Eq. (6) in the form

where $\hat{h}(\mathbf{r},t)=\hat{h}_{0}(\mathbf{r})+\hat{w}(\mathbf{r},t)$. In this equation, the explicit exchange term (the last term on the right) describes only the exchange with the ionic states that are not included in the basis, and the exchange with the basis states is hidden in the remaining terms due to the use of the CCC symmetry.

Here we introduce the matrix of inter-channel potentials

The first term on the right depends only on $\mathbf{r}_{1}$, the second - only on $t$.

We define the exchange operator for an arbitrary function as

For the function $\tilde{\psi}_{m}(\mathbf{r}_{2},t)$, due to its orthogonality to the ionic basis states, only terms with inter-electronic and external potentials remain in the exchange operator

The CCC symmetrization procedure can be described as the operator $\hat{C}$, the action of which on the function is expressed as

To represent the radial dependence of the wave function, discrete variable representation (DVR) was used based on the finite elements of the Gauss-Lobatto quadrature. The angular dependence was represented by a spherical harmonic expansion [6].

The time evolution in Eq. (11) was carried out using an implicit fourth-order scheme [7], which is a generalization of the Crank-Nicholson scheme. The system of linear equations to which this scheme leads was solved using the biconjugate gradient method with a preconditioner. Each multiplication by the Hamiltonian matrix was accompanied by the CCC symmetrization (15) before and after the multiplication, which ensured the CCC symmetry of the function.

Suppression of the unphysical reflection from the radial grid boundary was carried out using exterior complex scaling (ECS) for the radial variable.

Ionization amplitudes were extracted from the calculated wave function using the t-SURFFc method [8].

In the first test, we check the convergence of the energy of the ground states as the number of configurations $N_{s}$ in Eq. (4) increases. The energy of the ground states was calculated by evolution in imaginary time. As the two-electron targets, we consider the He and Li atoms as well as the H₂ molecule. The atoms in the present molecular formalism correspond to the limit of the inter-nuclear distance set to zero. Our test results in comparison with the literature values are shown in Table 1. It can be seen that adding more ionic states results in only a small change in the calculated ground state energies. This is a well-known effect for variational calculations: if the selected basis functions do not have the same set of singularities as the exact function, then the first few functions make the main contribution to the energy, but a lot of functions are needed to converge to the exact value. We can’t use very many ionic states in calculations without greatly slowing down the program. But we are willing to pay a fee in the form of inaccuracies in the energy of stationary states to ensure the correctness of the wave function after ionization.

In the second test, we calculate the photoelectron spectra driven by a short broad-band XUV pulse. These spectra, displayed in Fig. 1, show clearly the sequence of spN⁺, $N\geq 2$ autoionizing resonances when a multi-configuration expansion with $N_{s}=5$ is included in Eq. (4). From these spectra we obtained the resonance energies $E_{e}$, which are shown in Table 2.

We also calculated the lifetimes of the meta-stable auto-ionizing states corresponding to the resonances. To do this, we calculated the time dependence of the probability of finding an electron inside the simulation boundary ($r<r_{ECS}=63$ a.u.) after exposure to a long XUV pulse precisely tuned to resonance. Such a pulse causes the population of the corresponding doubly excited state. After the leakage of electrons emitted due to direct ionization, the probability decreases due to the decay of the meta-stable state, as clear from Fig. 2. So, using an exponential approximation of the time dependence, one can obtain the lifetimes of the meta-stable state. The calculated lifetimes are shown in Table 2 along with the state energies. They are in good agreement with the corresponding exact values shown in the same Table.