Theoretical methods for attosecond electron and nuclear dynamics: applications to the H₂ molecule

The total non-relativistic Hamiltonian for the internal motion of a diatomic molecule with two nuclei A and B, with masses M_A and M_B and nuclear charges Z_A and Z_B, and n electrons whose coordinates are referred to the geometrical center of the nuclei, can be expressed as

$$\begin{eqnarray}\hat{{\mathcal{H}}} & = & \displaystyle \sum _{i=1}^{n}\left[-\displaystyle \frac{1}{2}{{\rm{\nabla }}}_{i}^{2}-\displaystyle \frac{{Z}_{A}}{| {{\bf{r}}}_{i}-{{\bf{R}}}_{A}| }-\displaystyle \frac{{Z}_{B}}{| {{\bf{r}}}_{i}-{{\bf{R}}}_{B}| }\right]\\ & & +\displaystyle \frac{{Z}_{A}{Z}_{B}}{R}-\displaystyle \frac{1}{2\mu }{{\rm{\nabla }}}_{{\bf{R}}}^{2}-\displaystyle \sum _{i=1}^{n}\displaystyle \frac{1}{8\mu }{{\rm{\nabla }}}_{i}^{2}-\displaystyle \sum _{i\gt j}\displaystyle \frac{1}{2{M}_{T}}{{\rm{\nabla }}}_{i}{{\rm{\nabla }}}_{j}\\ & & -\displaystyle \frac{1}{2}\left(\displaystyle \frac{1}{{M}_{A}}-\displaystyle \frac{1}{{M}_{B}}\right)\displaystyle \sum _{i=1}^{n}{{\rm{\nabla }}}_{{\bf{R}}}{{\rm{\nabla }}}_{i},\end{eqnarray} \tag{ 6 }$$

where M_T = M_A+M_B is the total nuclear mass, μ = ${M}_{A}{M}_{B}/{M}_{T}$ is the reduced nuclear mass and R = $| {\bf{R}}|$ = $| {{\bf{R}}}_{A}-{{\bf{R}}}_{B}|$ is the internuclear distance. For a homonuclear molecule, e.g., H₂, M_A = M_B, the last term in equation (6) is strictly zero. As ${M}_{T},\mu \gg 1$, the mass polarization terms $\frac{1}{2{M}_{T}}{{\rm{\nabla }}}_{i}{{\rm{\nabla }}}_{j}$ and the $\frac{1}{8\mu }{{\rm{\nabla }}}_{i}^{2}$ terms can be neglected. Furthermore, we can now fix the molecular axis to reduce the size of the problem by removing the two Euler angles that describe the orientation of that axis in space. By further removing the three coordinates associated to the total centre of mass, the description of the quantum system remains with $3n+1$ coordinates, that is the internuclear distance R plus the $3n$ coordinates of the n electrons (denoted $\{{\bf{r}}\}$ for short). Henceforth we will deal with the following approximate Hamiltonian

$$\begin{eqnarray}&&{\hat{{\mathcal{H}}}}^{0}({\bf{r}};R)={\hat{T}}_{N}(R)+{\hat{{\mathcal{H}}}}_{{el}}({\bf{r}};R),\end{eqnarray} \tag{ 7 }$$

where ${\hat{T}}_{N}(R)$ = $-\frac{1}{2\mu }{{\rm{\nabla }}}_{R}^{2}$ is the relative kinetic energy of the nuclei and ${\hat{{\mathcal{H}}}}_{{el}}={\sum }_{i=1}^{n} [-{{\rm{\nabla }}}_{i}^{2}/2-{Z}_{A}/ | {{\bf{r}}}_{i}-{{\bf{R}}}_{A}| -{Z}_{B}/| {{\bf{r}}}_{i}-{{\bf{R}}}_{B}| ]\;+ ({Z}_{A}{Z}_{B})/R$ is the electronic Hamiltonian with a parametric dependence on the nuclear coordinate R.

The rotational degrees of freedom can be simply neglected when considering the timescales of interest in the sub-fs regime. Indeed, the rotational motion occurs in the picosecond regime, thus it is approximately three orders of magnitude slower than vibrations (typically taking place in tens of fs) and five to six orders of magnitude slower than electrons (attosecond regime). More importantly, the availability of kinematically complete experiments [99] along with ultrashort laser pulses of attosecond duration [6] has opened the way to work with fixed-in-space molecules. Indeed, in these experiments, the axial recoil approximation [100] is used to analyze the results since, after photo-absorption in the continuum, the time employed by the molecule to dissociate is much shorter than the rotational period of the molecule. In practice, it means that in dissociative ionization processes the set of vectors (${{\bf{K}}}_{{A}^{+}}$, ${{\bf{k}}}_{{e}^{-}}$, ${{\boldsymbol{\epsilon }}}_{p}$), where ${{\bf{K}}}_{{A}^{+}}$ is the momentum of the A⁺ ion, ${{\bf{k}}}_{{e}^{-}}$ is the electron momentum and ${{\boldsymbol{\epsilon }}}_{p}$ is the polarization direction, are measured in coincidence and can therefore be assigned to the body-fixed-frame triplet of vectors (${\bf{R}}$, ${{\bf{k}}}_{{e}^{-}}$, ${{\boldsymbol{\epsilon }}}_{p}$) for a given nuclear orientation of the molecular axis ${\bf{R}}$. To fully solve a quantum system with $3n+1$ independent coordinates is still a very demanding task, even for the seven-dimensional H₂ problem $({{\bf{r}}}_{1},{{\bf{r}}}_{2},R)$.

The wave function that describes the state of a molecule when interacting with a laser pulse is given by the solution of the TDSE:

$$\begin{eqnarray}&&i\displaystyle \frac{\partial }{\partial t}{\rm{\Psi }}(t)=\hat{{\mathcal{H}}}(t){\rm{\Psi }}(t),\end{eqnarray} \tag{ 8 }$$

where $\hat{{\mathcal{H}}}(t)$ is the total Hamiltonian, $\hat{{\mathcal{H}}}(t)$ = ${\hat{{\mathcal{H}}}}^{0}$+${\hat{V}}_{L}^{I}(t)$. Here ${\hat{{\mathcal{H}}}}^{0}$ represents the unperturbed Hamiltonian of the molecule described by equation (7), and ${\hat{V}}_{L}^{I}(t)$ is the time-dependent laser–molecule interaction potential. The TDSE must be solved by assuming that the molecular wave function is known at time t = t₀, i.e., before the laser–molecule interaction takes place. In this paper, the TDSE will be solved using the spectral method, in which the wave function ${\rm{\Psi }}(t)$ is expanded in the basis of eigenstates resulting from the diagonalization of the unperturbed Hamiltonian ${\hat{{\mathcal{H}}}}^{0}$ (the so-called close-coupling expansion). This approach has been widely used in atoms [101, 102] to describe ionization processes, even in the presence of autoionizing states. Usually, the integral over eigenstates representing the ionization continuum is replaced by a sum over discretized continuum states that result from the diagonalization of ${\hat{{\mathcal{H}}}}^{0}$ in a box. In general, the box must be large enough to provide a dense enough spectrum in the continuum that allows one to unambiguously resolve resonant states. In molecules, the practical implementation of this idea is not straightforward due to the presence of the nuclear degrees of freedom. Indeed, in general, electronic and nuclear degrees of freedom cannot be separated, unless in the framework of the BO approximation. However, this is not a good approximation when resonances are present, since the autoionization lifetime of these resonances can be comparable to the time needed by the nuclei to move. Thus, the general assumption behind the BO approximation, which is that electrons are much faster than nuclei, no longer holds in this context. Therefore, in order to obtain a meaningful basis of molecular eigenstates, one should diagonalize the whole molecular Hamiltonian ${\hat{{\mathcal{H}}}}^{0}$, not only the electronic Hamiltonian ${\hat{{\mathcal{H}}}}_{{el}}$. As a consequence, in the molecular case one must expand the time-dependent molecular wave function in a basis of vibronic states, in which electronic resonances are somewhat diluted and difficult to identify. Another difficulty arises when the autoionization decay time is very long, since in this case dissociation may occur well before the electron is ejected. In other words, a state that could potentially autoionize might not do it and behave as a truly bound state that dissociates into neutral species. This dual behavior of resonances does not exist in atoms because there are no sufficiently fast competing processes that can suppress autoionization.

To avoid diagonalizing the total molecular Hamiltonian ${\hat{{\mathcal{H}}}}^{0}$, we will adapt the stationary Feshbach theory for the solution of the TDSE. The original Feshbach projection method [98] was introduced in atoms by O'Malley et al (see [103] and references therein) and its extension to molecular systems was developed later [66, 67], following the seminal works of Bardsley [104] and Hazi, and Rescigno and Kurilla [105]. The Feshbach formalism is based on the introduction of projection operators ${\mathcal{P}}$ and ${\mathcal{Q}}$, satisfying completeness (${\mathcal{P}}$+${\mathcal{Q}}$ = 1), idempotency (${{\mathcal{P}}}^{2}$ = ${\mathcal{P}}$, ${{\mathcal{Q}}}^{2}$ = ${\mathcal{Q}}$) and orthogonality (${\mathcal{P}}{\mathcal{Q}}$ = ${\mathcal{Q}}{\mathcal{P}}$ = 0), which project the total molecular wave function onto non-resonant scattering-like and bound-like states, respectively. By using these projection operators, the molecular time-dependent wave function ${\rm{\Psi }}({\bf{r}},R,t)\;\equiv$ $\langle {\bf{r}},R| {\rm{\Psi }}(t)\rangle$ (in the combined space of electronic and nuclear coordinates) can be split off and represented by the state

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle ={\mathcal{Q}}| {\rm{\Psi }}(t)\rangle +{\mathcal{P}}| {\rm{\Psi }}(t)\rangle .\end{eqnarray} \tag{ 9 }$$

By inserting the latter equation in the TDSE (equation (8)) and projecting separately in both orthogonal subspaces, one obtains two time-dependent coupled equations, that in matrix form read

$$\begin{eqnarray}&&i\displaystyle \frac{\partial }{\partial t}\left[\begin{array}{c}{\mathcal{Q}}| {\rm{\Psi }}(t)\rangle \\ {\mathcal{P}}| {\rm{\Psi }}(t)\rangle \end{array}\right]=\left(\begin{array}{cc}{\mathcal{Q}}\hat{{\mathcal{H}}}{\mathcal{Q}} & {\mathcal{Q}}\hat{{\mathcal{H}}}{\mathcal{P}}\\ {\mathcal{P}}\hat{{\mathcal{H}}}{\mathcal{Q}} & {\mathcal{P}}\hat{{\mathcal{H}}}{\mathcal{P}}\end{array}\right)\left[\begin{array}{c}{\mathcal{Q}}| {\rm{\Psi }}(t)\rangle \\ {\mathcal{P}}| {\rm{\Psi }}(t)\rangle \end{array}\right].\end{eqnarray} \tag{ 10 }$$

One of the advantages of writing the TDSE in this way is that non-adiabatic effects are expected to be small within the resonant ${\mathcal{Q}}$ and non-resonant ${\mathcal{P}}$ subspaces. If one neglects non-adiabatic effects, the nuclear kinetic energy and the projection operators commute such that $[{\hat{T}}_{N}(R),{\mathcal{Q}}]$ = $[{\hat{T}}_{N}(R),{\mathcal{P}}]$ = 0 and ${\mathcal{Q}}{\hat{T}}_{N}(R){\mathcal{P}}$ = ${\mathcal{P}}{\hat{T}}_{N}(R){\mathcal{Q}}$ = 0. This means that the resonance–continuum coupling is entirely due to the electron–electron interaction, i.e., to ${\mathcal{Q}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{P}}$ (or its Hermitian conjugate ${\mathcal{P}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{Q}}$). This is usually a very good approximation since autoionization lifetimes of doubly excited states are almost exclusively due to the electron–electron interaction. Using these commutation rules, the total molecular Hamiltonian can be further separated in a block-diagonal unperturbed Hamiltonian [${\hat{{\mathcal{H}}}}_{{Fesh}}^{0}$ = ${\mathcal{Q}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{Q}}$ ⊕ ${\mathcal{P}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{P}}$] plus a total interaction potential [${\hat{{\mathcal{V}}}}_{L}^{I}(t)$ = ${\mathcal{Q}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{P}}$+${\mathcal{P}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{Q}}$+${\hat{V}}_{L}^{I}(t)$], which contains the explicit time-dependent term due to the laser–molecule interaction and two time-independent coupling terms that account for autoionization:

$$\begin{eqnarray}&&i\displaystyle \frac{\partial }{\partial t}\left[\begin{array}{l}{\mathcal{Q}}| {\rm{\Psi }}(t)\rangle \\ {\mathcal{P}}| {\rm{\Psi }}(t)\rangle \end{array}\right]\\ =\;\left\{\left(\begin{array}{cc}{\mathcal{Q}}[{\hat{T}}_{N}(R)+{\hat{{\mathcal{H}}}}_{{el}}]{\mathcal{Q}} & 0\\ 0 & {\mathcal{P}}[{\hat{T}}_{N}(R)+{\hat{{\mathcal{H}}}}_{{el}}]{\mathcal{P}}\end{array}\right)\right.\\ \quad +\left.\left(\begin{array}{cc}{\mathcal{Q}}{\hat{V}}_{L}^{I}(t){\mathcal{Q}} & {\mathcal{Q}}[{\hat{{\mathcal{H}}}}_{{el}}+{\hat{V}}_{L}^{I}(t)]{\mathcal{P}}\\ {\mathcal{P}}[{\hat{{\mathcal{H}}}}_{{el}}+{\hat{V}}_{L}^{I}(t)]{\mathcal{Q}} & {\mathcal{P}}{\hat{V}}_{L}^{I}(t){\mathcal{P}}\end{array}\right)\right\}\\ &&\quad \times \;\left[\begin{array}{c}{\mathcal{Q}}| {\rm{\Psi }}(t)\rangle \\ {\mathcal{P}}| {\rm{\Psi }}(t)\rangle \end{array}\right].\end{eqnarray} \tag{ 11 }$$

We use a spectral method to solve this coupled TDSE by expanding the total time-dependent wave function in terms of the eigenstates of the unperturbed Feshbach Hamiltonian ${\hat{{\mathcal{H}}}}_{{\rm{Fesh}}}^{0}$ = ${\mathcal{Q}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{Q}}$⊕${\mathcal{P}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{P}}$, from which we can build up a complete basis of vibronic states. Following common usage in theoretical chemistry, each vibronic state is written as a product of an electronic wave function ${\psi }_{n}({\bf{r}},R)$ describing the nth electronic state and a one-dimensional vibrational function ${\chi }_{{v}_{n}}(R)$, where v_n is the vibrational quantum number associated to the nth potential energy curve. We note that, for dissociative states, v_n is a continuous index. Methods to obtain the vibrational wave functions are discussed in detail in section 3.5. The complete set of vibronic states includes (i) the BO wave functions $\{{\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}({\bf{r}},R)\cdot {\chi }_{{v}_{\alpha }}(R)\}$ with total energies $\{{W}_{\alpha {v}_{\alpha }}={\varepsilon }_{\alpha }+{E}_{{v}_{\alpha }}\}$, which are approximate solutions of the ${\mathcal{P}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{P}}$ eigensystem (α denotes a given ionic threshold, ${\varepsilon }_{\alpha }$ is the electron kinetic energy in excess above the ionic threshold α, ${{\ell }}_{\alpha }$ is its angular momentum and ${E}_{{v}_{\alpha }}$ is the vibrational energy of state ${v}_{\alpha }$ in the potential energy curve of the electronic state α); (ii) the BO wave functions $\{{\psi }_{r}({\bf{r}},R)\cdot {\chi }_{{v}_{r}}(R)\}$ with total energies ${W}_{{{rv}}_{r}}$, which are approximate solutions of the ${\mathcal{Q}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{Q}}$ eigensystem; and (iii) the vibronic states associated to the electronic bound states $\{{\psi }_{b}({\bf{r}},R)\cdot {\chi }_{{v}_{b}}(R)\}$ with total energies $\{{W}_{{{bv}}_{b}}\}$. These bound states mostly pertain to the ${\mathcal{P}}$ subspace. However, in practice, they are more conveniently computed without the Feshbach partitioning.

The expansion of the total time-dependent wave function in our basis of BO vibronic states is written in the Schrödinger picture as:

$$\begin{eqnarray}{\rm{\Psi }}({\bf{r}},R,t) & = & \displaystyle \sum _{b}{\displaystyle \sumint }_{{v}_{b}}{C}_{{{bv}}_{b}}(t){\psi }_{b}({\bf{r}},R){\chi }_{{v}_{b}}(R){e}^{-{{iW}}_{{{bv}}_{b}}t}\\ & & +\displaystyle \sum _{r}{\displaystyle \sumint }_{{v}_{r}}{C}_{{{rv}}_{r}}(t){\psi }_{r}({\bf{r}},R){\chi }_{{v}_{r}}(R){e}^{-{{iW}}_{{{rv}}_{r}}t}\\ & & +\displaystyle \sum _{\alpha {{\ell }}_{\alpha }}\displaystyle \int {\rm{d}}{\varepsilon }_{\alpha }{\displaystyle \sumint }_{{v}_{\alpha }}{C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}(t){\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}({\bf{r}},R)\\ & & \times {\chi }_{{v}_{\alpha }}(R){e}^{-{{iW}}_{\alpha {v}_{\alpha }}t}.\end{eqnarray} \tag{ 12 }$$

To solve the TDSE defined in equation (8), we introduce the latter ansatz in equation (11) to obtain a system of coupled linear differential equations that can be written in compact matrix-vector form as (in the interaction picture)

$$\begin{eqnarray}&&i\displaystyle \frac{\partial }{\partial t}\left(\begin{array}{c}{{\bf{C}}}_{{{bv}}_{b}}^{I}\\ {{\bf{C}}}_{{{rv}}_{r}}^{I}\\ {{\bf{C}}}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{I;{\varepsilon }_{\alpha }}\end{array}\right)\\ =\;\left(\begin{array}{ccc}{{\bf{V}}}_{L}^{I}{(t)}_{{{bv}}_{b}}^{b^{\prime} {v}_{{b}^{\prime }}^{\prime }} & {{\bf{V}}}_{L}^{I}{(t)}_{{{bv}}_{b}}^{{{rv}}_{r}} & {{\bf{V}}}_{L}^{I}{(t)}_{{{bv}}_{b}}^{\alpha {{\ell }}_{\alpha }{v}_{\alpha }{\varepsilon }_{\alpha }}\\ {{\bf{V}}}_{L}^{I}{(t)}_{{{rv}}_{r}}^{{{bv}}_{b}} & {{\bf{V}}}_{L}^{I}{(t)}_{{{rv}}_{r}}^{r^{\prime} {v}_{{r}^{\prime }}^{\prime }} & \begin{array}{l}\left[{\mathcal{P}}{{\mathcal{H}}}_{{el}}{\mathcal{Q}}\right.\\ +\left.{{\bf{V}}}_{L}^{I}(t)\right]{\ }_{{{rv}}_{r}}^{\alpha {{\ell }}_{\alpha }{v}_{\alpha }{\varepsilon }_{\alpha }}\end{array}\\ {{\bf{V}}}_{L}^{I}{(t)}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }{\varepsilon }_{\alpha }}^{{{bv}}_{b}} & \begin{array}{l}\left[{\mathcal{Q}}{{\mathcal{H}}}_{{el}}{\mathcal{P}}\right.\\ +\left.{{\bf{V}}}_{L}^{I}(t)\right]{\ }_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }{\varepsilon }_{\alpha }}^{{{rv}}_{r}}\end{array} & {{\bf{V}}}_{L}^{I}{(t)}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }{\varepsilon }_{\alpha }}^{{\alpha }^{\prime }{{\ell }}_{{\alpha }^{\prime }}^{\prime }{v}_{{\alpha }^{\prime }}^{\prime }{\varepsilon }_{{\alpha }^{\prime }}^{\prime }}\end{array}\right)\\ &&\quad \times \;\left(\begin{array}{c}{{\bf{C}}}_{{{bv}}_{b}}^{I}\\ {{\bf{C}}}_{{{rv}}_{r}}^{I}\\ {{\bf{C}}}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{I;{\varepsilon }_{\alpha }}\end{array}\right).\end{eqnarray} \tag{ 13 }$$

The couplings due to the time-dependent laser interaction, ${{\bf{V}}}_{L}^{I}(t)$, correspond to the dipole matrix elements in the long-wave approximation ($e{\bf{r}}$ in length gauge and $\frac{e}{{mc}}{\bf{p}}$ in velocity gauge) times the dependence on the laser pulse, i.e. ${{\bf{V}}}_{L}^{I}(t)$ = $e{\bf{r}}\cdot {\bf{E}}(t)$ in the length gauge and ${{\bf{V}}}_{L}^{I}(t)$ = $\frac{e}{{mc}}{\bf{p}}\cdot {\bf{A}}(t)$ in the velocity gauge. The crossed terms ${\mathcal{Q}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{P}}$+${\mathcal{P}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{Q}}$ are active during and after the interaction with the laser pulse. For atoms, this term becomes zero when the resonance (doubly excited) states are fully depleted, so that ${\mathrm{lim}}_{t\to \infty }{\mathcal{Q}}{\rm{\Psi }}(t)$ = 0 and ${\mathrm{lim}}_{t\to \infty }{\mathcal{P}}{\rm{\Psi }}(t)$ = ${\rm{\Psi }}(t)$ (for an application of the present TDFCC method to atoms see [106]). However, this is no longer true for molecules, because the time needed to dissociate into neutrals may be shorter than the autoionization lifetimes, and as a consequence part of the wave function norm can remain in the ${\mathcal{Q}}$ subspace. In other words, ${\mathrm{lim}}_{t\to \infty }{\mathcal{Q}}{\rm{\Psi }}(t)$ = ${{\rm{\Phi }}}_{{diss}}$, where ${{\rm{\Phi }}}_{{diss}}$ is the asymptotic representation of the dissociative channels in which no ionization has occurred. As in atomic systems, ${\mathrm{lim}}_{t\to \infty }{\mathcal{P}}{\rm{\Psi }}(t)$ = ${{\rm{\Phi }}}_{{ion}}$, where ${{\rm{\Phi }}}_{{ion}}$ is the asymptotic representation of the ionization channels. Asymptotically,

$$\begin{eqnarray}&&\underset{t\to \infty }{\mathrm{lim}}{\hat{{\mathcal{H}}}}^{0}={\mathcal{Q}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{Q}}+{\mathcal{P}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{P}},\end{eqnarray} \tag{ 14 }$$

so that, in that limit, the eigenfunctions of ${\mathcal{Q}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{Q}}$ and ${\mathcal{P}}{\hat{{\mathcal{H}}}}^{0}{\mathcal{P}}$ are also eigenfunctions of the total Hamiltonian. Therefore, the expansion coefficients ${C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}$ of the ${\mathcal{P}}$ subspace truly represent physical ionization amplitudes once the system reaches its stationary asymptotic limit. The most attractive feature of the present time-dependent Feshbach formulation is that it provides a straightforward interpretation of autoionization processes in time, in contrast to the Feshbach stationary approach [80]. In the time-dependent formulation, the total scattering wave function is built up during time propagation, which allows for the simultaneous excitation of the resonant ${\mathcal{Q}}| {\rm{\Psi }}\rangle$ and non-resonant ${\mathcal{P}}| {\rm{\Psi }}\rangle$ subspaces from the initial ground state and their mixing through the interaction potential ${\hat{{\mathcal{V}}}}_{L}^{I}(t)$. Again, since we are using a representation of the time-dependent wave function in a basis of molecular states, the transition amplitudes are directly related to the expansion coefficients in equation (12) at $t=\infty$.

The present spectral method has an added advantage regarding computational effort, since the calculation of dipole and ${\mathcal{Q}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{P}}$ coupling matrix elements included in the interaction potential terms ${\hat{{\mathcal{V}}}}_{L}^{I}(t)$ in equation (13) is only performed once. The specific laser parameters entering in ${\hat{{\mathcal{V}}}}_{L}^{I}(t)$ appear as trivial factors, so that in practice one must only invest computer effort in solving the system of coupled linear differential equations (13) for each specific laser field. This is very convenient, because the calculation of the coupling matrix elements is by far the most demanding part from the computational point of view. Indeed, the evaluation of these couplings involves electronic and vibrational states and integrals between them, accounting for both the ${\bf{r}}$ and R dependences. For instance, the matrix elements in the last column of the matrix representation of the TDSE in equation (13) represent the coupling between final vibronic states, $\{{\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}\cdot {\chi }_{{v}_{\alpha }}\}$ and (i) initial or intermediate vibronic bound states $\{{\psi }_{b}\cdot {\chi }_{{v}_{b}}\}$ (upper block), (ii) intermediate vibronic resonant states $\{{\psi }_{r}\cdot {\chi }_{{v}_{r}}\}$ (middle block) and (iii) intermediate electronic continuum states $\{{\psi }_{{\alpha }^{\prime }{{\ell }}_{{\alpha }^{\prime }}^{\prime }}^{{\varepsilon }_{{\alpha }^{\prime }}^{\prime }}\cdot {\chi }_{{v}_{{\alpha }^{\prime }}^{\prime }}\}$ that account for above threshold ionization processes (lower block). These couplings have the general form

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{{R}_{\mathrm{max}}}{\rm{d}}{R}{\chi }_{{v}_{n}}(R)\cdot \langle {\psi }_{n}| {\hat{{\mathcal{V}}}}_{L}^{I}(t)| {\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}\rangle (R)\cdot {\chi }_{{v}_{\alpha }}(R),\end{eqnarray} \tag{ 15 }$$

where n represents the labels required for the states in any of the above-mentioned cases and $[0,{R}_{\mathrm{max}}]$ indicates the extension of the nuclear radial box. To compute these vibronic couplings, it is then required to calculate first the electronic matrix elements $\langle {\psi }_{n}| {\hat{{\mathcal{V}}}}_{L}^{I}(t)| {\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}\rangle (R)$ in a dense grid of internuclear distances within the interval $[0,{R}_{\mathrm{max}}]$, namely, the dipolar transition matrix elements $\langle {\psi }_{n}| {\hat{\mu }}_{L,V}| {\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}\rangle$ with ${\hat{\mu }}_{L}$ = $e{\bf{r}}$ (length) or ${\hat{\mu }}_{V}$ = $\frac{e}{{mc}}{\bf{p}}$ (velocity) as well as the electrostatic Feshbach couplings $\langle {\psi }_{n}| {\mathcal{Q}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{P}}| {\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}\rangle$ responsible for the resonance decay into the continuum. Consequently, suitable electronic and nuclear structure calculations are needed to generate the set of vibronic states and their energies $\{{\psi }_{n},{\chi }_{{v}_{n}},{W}_{{{nv}}_{n}}\}$. The calculation of vibronic states for electronic bound states is rather standard, while it is not so trivial for continuum states, especially when both the electronic and the nuclear continuum are involved.

In the following, we describe the method that we have used to obtain the vibronic states used in the close-coupling expansion of the time-dependent wave function. Readers not interested in these technical details may skip this section and go directly to section 4, where a number of different applications of the TDFCC method are discussed.

For the computation of n-electron (bound, resonant and continuum) molecular states of a diatomic molecule, one can make use of a configuration interaction (CI) method based on expansions in terms of antisymmetrized products of one-particle wave functions. For this, the molecular orbitals for the one-electron molecular ion (${Z}_{A},{e}^{-},{Z}_{B}$) must firstly be computed. Consistently with equation (6), they are eigensolutions of the Schrödinger equation

$$\begin{eqnarray}&&\left[-\displaystyle \frac{1}{2}{{\boldsymbol{\nabla }}}_{{\bf{r}}}^{2}-\displaystyle \frac{{Z}_{A}}{| {\bf{r}}-{{\bf{R}}}_{A}| }-\displaystyle \frac{{Z}_{B}}{| {\bf{r}}-{{\bf{R}}}_{B}| }-{\epsilon }_{n}(R)\right]{\varphi }_{n}({\bf{r}};R)=0.\end{eqnarray} \tag{ 16 }$$

This molecular two-center Coulomb system is a well known exactly solvable problem in molecular physics. It represents, e.g., the ${{\rm{H}}}_{2}^{+}$ molecule (Z_A = Z_B = 1). The Schrödinger equation (16) is fully separable in prolate spheroidal coordinates (see for instance [107]) and both energies and wave functions can be computed with arbitrary precision [108–110]. This one-electron molecular Hamiltonian, which we call ${\hat{h}}_{{el}}$, remains invariant under the rotation of the internuclear axis for both heteronuclear (point group ${C}_{\infty v}$) and homonuclear molecules (point group ${D}_{\infty h}$) and, consequently, it commutes with the operator ${\hat{l}}_{z}$ such that the molecular orbitals can simultaneously be eigenstates of ${\hat{h}}_{{el}}$ and ${\hat{l}}_{z}$. The value of the quantum number λ = $| m|$ corresponds to the absolute value of the eigenvalue of the angular momentum operator ${\hat{l}}_{z}$ along the internuclear axis Z.

The use of the exact one-electron diatomic molecular (OEDM) orbitals in an n-electron configuration interaction approach along with a scheme well suited to describe the electronic continuum is rather cumbersome [97]. To avoid the difficulties inherent to using non-${{\mathcal{L}}}^{2}$ states and to recover the essence of the partial wave analysis common to many scattering problems, the one-electron wave functions can be expressed in terms of a one-center spherical expansion. For a fixed value of m, the expansion reads

$$\begin{eqnarray}&&{\varphi }_{n}^{\lambda ,m}({\bf{r}};R)=\displaystyle \sum _{{\ell }=0}^{{{\ell }}_{\mathrm{max}}}\displaystyle \frac{{U}_{n{\ell }}(r;R)}{r}{{\mathcal{Y}}}_{{\ell }}^{m}(\theta ,\phi ),\end{eqnarray} \tag{ 17 }$$

where ${{\mathcal{Y}}}_{{\ell }}^{m}$ refers to the spherical harmonics and we have kept the value of m with its sign to make clear that different values of m must be combined with their respective signs when preparing n-electron configurations. The radial wave function ${U}_{n{\ell }}(r;R)$ can be expanded in terms of suitable basis functions for each internuclear distance R. Our choice for the basis set is B-splines and, therefore, ${U}_{n{\ell }}(r;R)$ = ${\displaystyle \sum }_{i=1}^{{N}_{B}}{c}_{i{\ell }}^{n}(R){B}_{i}^{k}(r)$, where N_B is the number of B-spline functions. B-splines are piecewise polynomials of order k defined inside a radial box of length ${r}_{\mathrm{max}}$ [111]. The set of k B-splines at the electron position r is constructed by defining a knot sequence ${\{{t}_{i}\}}_{i=1}^{{N}_{B}+k}$ within the electronic radial box with $r\in [0,{r}_{\mathrm{max}}]$, using the recursion formula [111, 112]

$$\begin{eqnarray}&&{B}_{i}^{k}(r)=\displaystyle \frac{r-{t}_{i}}{{t}_{i+k-1}-{t}_{i}}{B}_{i}^{k-1}(r)+\displaystyle \frac{{t}_{i+k}-r}{{t}_{i+k}-{t}_{i+1}}{B}_{i+1}^{k-1}(r);\end{eqnarray} \tag{ 18 }$$

and setting the B-splines of order k = 1 as ${B}_{i}^{1}(r)=1$ for ${t}_{i}\leqslant r\lt {t}_{i+1}$, and ${B}_{i}^{1}(r)=0$ otherwise. The two-boundary eigenvalue problem is solved subject to the boundary conditions ${U}_{n{\ell }}(0)$ = ${U}_{n{\ell }}({r}_{\mathrm{max}})$ = 0, which are automatically satisfied by removing the first ${B}_{1}^{k}(r)$ and the last ${B}_{{N}_{B}}^{k}(r)$ B-spline from the original set.

The latter expansion turns the solution of the original set of differential equations into a generalized eigenvalue problem $({\bf{h}}-{\boldsymbol{\epsilon }}{\bf{s}}){\bf{c}}$ = 0 for each internuclear distance R, where ${\bf{c}}$ is the vector for the expansion coefficients ${c}_{i{\ell }}^{n}(R)$. The matrix elements of the Hamiltonian and the overlap in terms of the basis ${\{{B}_{i}^{k}\cdot {{\mathcal{Y}}}_{{\ell }}^{m}\}}_{i=1,{\ell }=0}^{{N}_{B},{{\ell }}_{\mathrm{max}}}$, for fixed m and R, are

$$\begin{eqnarray}{{\bf{h}}}_{i{\ell },i^{\prime} {\ell }^{\prime} } & = & {\displaystyle \int }_{0}^{{r}_{\mathrm{max}}}{\rm{d}}{r}{\displaystyle \int }_{0}^{\pi }{\rm{d}}\theta \mathrm{sin}\theta {\displaystyle \int }_{0}^{2\pi }{\rm{d}}\phi {B}_{i}^{k}(r)\\ & & \times {{\mathcal{Y}}}_{{\ell }}^{m*}(\theta ,\phi ){\hat{h}}_{{el}}{B}_{{i}^{\prime }}^{k}(r){{\mathcal{Y}}}_{{\ell }^{\prime} }^{m}(\theta ,\phi )\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\bf{s}}}_{i{\ell },i^{\prime} {\ell }^{\prime} }={\delta }_{{\ell },{\ell }^{\prime} }{\displaystyle \int }_{0}^{{r}_{\mathrm{max}}}{B}_{i}^{k}(r){B}_{{i}^{\prime }}^{k}(r){\rm{d}}{r}\end{eqnarray} \tag{ 20 }$$

respectively. The size of the one-electron eigensystem is ${N}_{B}\times {{\ell }}_{\mathrm{max}}$.

In contrast to atoms, the use of a one-center wave function expansion to describe a molecular two-center problem imposes that the electron-nucleus Coulomb potential ${V}^{q}={-Z}_{q}/| {\bf{r}}-{{\bf{R}}}_{q}|$ is more conveniently expanded in spherical coordinates. When the origin of the one center expansion coincides with the nuclear center of mass, the resulting effective Hamiltonian expressed in spherical polar coordinates is

$$\begin{eqnarray}{\hat{h}}_{{el}} & = & -\displaystyle \frac{{d}^{2}}{{{dr}}^{2}}+\displaystyle \frac{{\hat{L}}^{2}}{2{r}^{2}}-\displaystyle \sum _{q=A}^{B}{Z}_{q}\displaystyle \sum _{{{\ell }}_{q}=0}^{\infty }\displaystyle \frac{4\pi }{2{{\ell }}_{q}+1}\displaystyle \frac{{r}_{\lt }^{{{\ell }}_{q}}}{{r}_{\gt }^{{{\ell }}_{q}+1}}\\ & & \times \displaystyle \sum _{{m}_{q}=-{{\ell }}_{q}}^{+{{\ell }}_{q}}{{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}*}({\theta }_{q},{\phi }_{q}){{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}}(\theta ,\phi ),\end{eqnarray} \tag{ 21 }$$

where ${r}_{\lt }$ = min($r,{R}_{q}$), ${r}_{\gt }$ = max($r,{R}_{q}$) and $({\theta }_{q},{\phi }_{q})$ are the spherical angles of nucleus q=$\{A,B\}$. Accordingly, the corresponding matrix elements for the one-electron Hamiltonian adopt the form

$$\begin{eqnarray}&&{{\bf{h}}}_{i,{\ell },i^{\prime} {\ell }^{\prime} }={{\bf{t}}}_{i{\ell },i^{\prime} {\ell }^{\prime} }+\displaystyle \sum _{q=\{A,B\}}{{\bf{v}}}_{i{\ell },i^{\prime} {\ell }^{\prime} }^{q},\end{eqnarray} \tag{ 22 }$$

where the kinetic energy term

$$\begin{eqnarray}{{\bf{t}}}_{i{\ell },i^{\prime} {\ell }^{\prime} } & = & \left(-\displaystyle \frac{1}{2}{\displaystyle \int }_{0}^{{r}_{\mathrm{max}}}{B}_{i}^{k}(r)\displaystyle \frac{{{\rm{d}}}^{2}}{{{\rm{d}}{r}}^{2}}{B}_{{i}^{\prime }}^{k}(r){\rm{d}}{r}\right.\\ & & +\left.\displaystyle \frac{{\ell }({\ell }+1)}{2}{\displaystyle \int }_{0}^{{r}_{\mathrm{max}}}\displaystyle \frac{{B}_{i}^{k}(r){B}_{{i}^{\prime }}^{k}(r)}{{r}^{2}}{\rm{d}}{r}\right){\delta }_{{\ell }{\ell }^{\prime} }\end{eqnarray} \tag{ 23 }$$

is block diagonal in ${\ell }{\rm{s}}$ and the matrix elements for the Coulomb interaction read

$$\begin{eqnarray}{{\bf{v}}}_{i{\ell },i^{\prime} {\ell }^{\prime} }^{q} & = & -{Z}_{q}\displaystyle \sum _{{{\ell }}_{q}=0}^{\infty }\displaystyle \frac{4\pi }{2{{\ell }}_{q}+1}{\displaystyle \int }_{0}^{{r}_{\mathrm{max}}}{B}_{i}^{k}(r)\displaystyle \frac{{r}_{\lt }^{{{\ell }}_{q}}}{{r}_{\gt }^{{{\ell }}_{q}+1}}{B}_{{i}^{\prime }}^{k}(r){\rm{d}}{r}\\ & & \times \displaystyle \sum _{{m}_{q}=-{{\ell }}_{q}}^{+{{\ell }}_{q}}{{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}}({\theta }_{q},{\phi }_{q})\\ & & \times {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{\pi }{{\mathcal{Y}}}_{{\ell }}^{m*}(\theta ,\phi ){{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}}(\theta ,\phi ){{\mathcal{Y}}}_{{\ell }^{\prime} }^{{m}^{\prime }}\;(\theta ,\phi )\mathrm{sin}\theta {\rm{d}}\theta {\rm{d}}\phi .\end{eqnarray} \tag{ 24 }$$

In the latter expressions all angular integrals are analytical (Gaunt formula) and the radial ones can be readily computed to machine accuracy using a Gauss–Legendre quadrature since the integrands are reduced to polynomials (except the Coulomb potential terms). Incidentally, since we adopt a geometry where the linear molecule is oriented along the Z-axis and the origin is located at the nuclear geometrical center in homonuclear molecules, the spherical harmonics for the nuclei ${{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}}({\theta }_{q},{\phi }_{q})$ adopt simple constant values for $({\theta }_{A},{\phi }_{A})$ = $(0,0)$ and $({\theta }_{B},{\phi }_{B})$ = $(\pi ,0)$, namely, ${{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}}({\theta }_{A}=0,{\phi }_{A})=0$ = ${(-1)}^{{{\ell }}_{q}}$ ${{\mathcal{Y}}}_{{{\ell }}_{q}}^{{m}_{q}}$ $({\theta }_{B}=\pi ,{\phi }_{B}=0)$ =$\sqrt{(2{{\ell }}_{q}+1)/4\pi }$ if m_q=0 and they are identically zero if ${m}_{q}\ne 0$. Given that the pair $\{{{\ell }}_{q},{m}_{q}\}$ satisfies $| {\ell }-{\ell }^{\prime} | \leqslant {{\ell }}_{q}\leqslant {\ell }+{\ell }^{\prime}$ and $m+m^{\prime} +{m}_{q}$ = 0 (with the only contribution from m_q = 0 and m = −m') the quality in the representation of the molecular Coulomb interaction depends on the number of angular momenta ${{\ell }}_{\mathrm{max}}$ introduced in expansion (17). Notably, this number must increase as the internuclear distance increases since departure from the atom-like spherical symmetry becomes more pronounced.

In the case of homonuclear molecules, Z_A = Z_B, the electronic wave function has an inversion symmetry, $\hat{i}[\varphi ({\bf{r}})]$ = $\varphi (-{\bf{r}})$ = $\pi \varphi ({\bf{r}})$, and it remains unchanged (π = +1 or gerade (g)) or it changes its sign (π = −1 or ungerade (u)) against inversion of the electronic coordinates. In our single center expansion, the inversion symmetry operation only affects the spherical harmonics, i.e., ${{\mathcal{Y}}}_{{\ell }}^{m}(-\hat{{\bf{r}}})$ = ${(-1)}^{{\ell }}{{\mathcal{Y}}}_{{\ell }}^{m}(\hat{{\bf{r}}})$ and, consequently, only even values of the angular momentum contribute to orbitals of gerade symmetry and odd values to orbitals of ungerade symmetry. Furthermore, because both the angular momentum operator ${\hat{l}}_{z}$ and the inversion parity operator $\hat{i}$ commute with the electronic Hamiltonian ${\hat{h}}_{{el}}$, separate calculations can be performed for the different values of λ = 0 (σ), λ = 1 (π), λ = 2 (δ), λ = 3 (ϕ), etc, and for gerade and ungerade symmetries, using angular momentum expansions with $\lambda \leqslant {\ell }\leqslant {{\ell }}_{\mathrm{max}}$. Typically, for the computation of ${\sigma }_{g}$ orbitals, one can use a set of 180 B-splines of order k=8 in a linear sequence of knot-points enclosed in an electronic radial box of length ${r}_{\mathrm{max}}$=60 a.u. with an expansion using only even partial waves in equation (17) from ${\ell }$=0 up to ${{\ell }}_{\mathrm{max}}$ = 16. Note that different angular momenta will contribute for each symmetry. To compute ${\pi }_{u}$ orbitals ($\lambda =1$) the partial waves included are ${\ell }=1,3,\ldots ,15$, while for ${\delta }_{u}$ molecular orbitals (λ = 2) the wave function is expanded in terms of seven odd partial waves with ${\ell }$ = 3, 5, ..., 15.

It is worth noting that, in spite of the advantages of using B-splines in the representation of the continuum [112], there are a few limitations on the reliability of the subsquent CI procedure. In fact, these orbitals fail to accurately reproduce the inner electronic density close to the nuclei around the equilibrium distance. To improve the description of bound states in the case of H₂, the CI representation of the inner wave function is usually complemented with additional Slater-type orbitals (STOs), i.e., functions of the type ${r}^{n}\mathrm{exp}(-{\gamma }_{n{\ell }}r)\cdot {{\mathcal{Y}}}_{{\ell }}^{m}(\hat{{\bf{r}}})$ (typically using n = 1–10, ${\ell }\;=$ 0–11 and ${\gamma }_{n{\ell }}$ = 2.8 for all n and ${\ell }$). For convenience, these supplementary STO basis functions are expanded in the same basis of B-splines. Despite the loss of strict orthogonality and the inclusion of the corresponding overlaps, this addition significantly improves the energy convergence for the ground state in H₂ close to the internuclear equilibrium distance R_e ∼1.4 a.u.

In principle, with the solution of the one-particle Hamiltonian we have at our disposal a full set of approximate OEDM orbitals $\{{\varphi }_{n}^{\lambda ,m}\}(R)$ and also their corresponding electronic energies $\{{\epsilon }_{n}^{\lambda }\}(R)$. In the case of H₂, the latter energies represent the ionization energies at different values of R. A clear disadvantage of one-center expansions is their slower convergence as one moves away from the spherical symmetry at R = 0; the eigenvalues can reproduce the exact OEDM energies [109] at large internuclear distances only at the cost of huge angular momentum expansions, which become impractical. This issue can still be circumvented. First, in most applications, e.g. in H₂, the molecule is initially in its ground state, with a density localized close to the equilibrium internuclear distance, and the computed dipole matrix elements can easily be converged up to a rather large value of R. Secondly, although the slow convergence at large internuclear distances may affect the R-asymptotic behavior and the energies of the vibronic wave functions, and by extension the continuum (dissociating) nuclear wave functions, the problem can be overcome by using accurate potential energy curves resulting from standard quantum chemistry calculations for the ground and singly excited states (see, e.g., [113] for the H₂ molecule) and/or the exact OEDM potential energy curves for the ionic thresholds.

The molecular non-resonant electronic continua are the eigenstates of the ${\mathcal{P}}$ projected electronic Hamiltonian. One has to solve the eigenvalue problem

$$\begin{eqnarray}&&[{\mathcal{P}}{\hat{{\mathcal{H}}}}_{{el}}{\mathcal{P}}-E]{\psi }_{\mu }^{{\varepsilon }_{\mu }}({\bf{r}};R)=0\end{eqnarray} \tag{ 25 }$$

for each internuclear distance R. Henceforth we will use a more compact notation and define a collective channel index, μ = (α, ${{\ell }}_{\alpha }$), encompassing the previously defined ionic threshold label α, with energy ${E}_{\alpha }(R)$, and the angular momentum of the scattered electron ${{\ell }}_{\alpha }$ (${{\ell }}_{\mu }$ from now on) that is compatible with the molecular symmetries of the $(n-1)$-electron target. Also, ${\varepsilon }_{\mu }\equiv {\varepsilon }_{\alpha }$.

The total electronic energy of the n-electron system satisfies E = ${E}_{\alpha }(R)$+${\varepsilon }_{\mu }$, where ${E}_{\alpha }(R)$ is the potential energy curve of the ionic threshold α (in the case of H₂, it is simply one of the OEDM energies from the set $\{{\epsilon }_{n}^{\lambda }\}(R)$ previously described) and ${\varepsilon }_{\mu }$ is the kinetic energy for the emitted electron in channel μ. In our implementation, instead of explicitly calculating the eigenstates of the projected Hamiltonian in equation (25), we use multichannel scattering theory to evaluate electronic continuum states from a basis that spans the ${\mathcal{P}}$ subspace. Thus we solve instead

$$\begin{eqnarray}&&\;[{\hat{{\mathcal{H}}}}_{{el}}-E]{\mathcal{P}}{{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}({\bf{r}},R)=0,\end{eqnarray} \tag{ 26 }$$

where ${\mathcal{P}}{{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}$ is the incoming scattering wave function that represents one electron being ionized from the molecule. It is constructed to be an eigenfunction of the set of operators $\hat{{\rm{\Omega }}}\equiv \{{\hat{S}}^{2},{\hat{S}}_{z},{\hat{L}}_{z},\hat{\sigma }\}$ (plus $\hat{i}$ in the case of a homonuclear molecule) with their corresponding eigenvalues ${\rm{\Omega }}\equiv \{S(S+1),{M}_{s},\pm {\rm{\Lambda }},\sigma \}$ (and π). S is the total electronic spin, M_s is its z-component, Λ the absolute value of the z-component of the total electronic angular momentum (${\rm{\Sigma }},{\rm{\Pi }},{\rm{\Delta }},...$), σ is the reflexion symmetry with respect to the XZ plane (±1) and π is the symmetry with respect to the inversion center (g for gerade and u for ungerade).

This projected wave function, for the channel μ and total electronic energy E, is built through a multichannel close-coupling expansion in the form (for notational simplicity, in the rest of the formulation we will consider the case of a two-electron molecule, so that spin functions can be factored out)

$$\begin{eqnarray}&&{\mathcal{P}}{{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}({{\bf{r}}}_{1},{{\bf{r}}}_{2},R)={\rm{\Theta }}\displaystyle \sum _{\mu ^{\prime} }^{{N}_{c}}\left[{\varphi }_{\mu ^{\prime} }^{{\rm{\Omega }}}({{\bf{r}}}_{1},{\hat{r}}_{2}){F}_{\mu {\mu }^{\prime }}^{{\rm{\Omega }}-}({r}_{2})\right],\end{eqnarray} \tag{ 27 }$$

where Θ is the symmetrization (antisymmetrization) operator for singlet (triplet) states, N_c is the number of open channels, and the channel functions, ${\varphi }_{{\mu }^{\prime }}^{{\rm{\Omega }}}$, correspond to target states of the hydrogen molecular ion coupled with the angular functions of the scattered electron such that they are eigenfuntions of the operator $\hat{{\rm{\Omega }}}$. In the more general case of an n-electron diatomic molecule, ${{\bf{r}}}_{1}$ represents the coordinates of the $n-1$ electrons that remain bound, and ${\varphi }_{{\mu }^{\prime }}^{{\rm{\Omega }}}({{\bf{r}}}_{1},{\hat{r}}_{2})$ includes the spin functions for all n electrons (also, Θ can only be the antisymmetrization operator).

Open and closed channels are dictated by the total electronic energy E; for open channels ${\varepsilon }_{\mu }$ = $E-{E}_{\alpha }(R)\geqslant 0$ and for closed channels ${\varepsilon }_{\mu }\leqslant 0$. Convergence of the close-coupling expansion can only be obtained by including both open and closed channels.

In particular, the radial wave function of the ionized electron ${F}_{\mu {\mu }^{\prime }}^{{\rm{\Omega }}-}(r)$ has the following (energy-normalized) asymptotic behaviour:

$$\begin{eqnarray}&&\underset{r\to \infty }{\mathrm{lim}}{F}_{\mu {\mu }^{\prime }}^{{\rm{\Omega }}-}(r)=\sqrt{\displaystyle \frac{1}{2\pi {k}_{\mu }}}\displaystyle \frac{1}{{ir}}\left({\delta }_{\mu {\mu }^{\prime }}{e}^{i{\theta }_{\mu }}-{{\mathcal{S}}}_{\mu {\mu }^{\prime }}^{{\rm{\Omega }}\dagger }{e}^{-i{\theta }_{{\mu }^{\prime }}}\right),\end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray}&&{\theta }_{\mu }={k}_{\mu }r-\displaystyle \frac{1}{2}{{\ell }}_{\mu }\pi +\displaystyle \frac{Z-1}{{k}_{\mu }}\mathrm{log}(2{k}_{\mu }r)+{\sigma }_{{{\ell }}_{\mu }}({k}_{\mu }),\end{eqnarray} \tag{ 29 }$$

where ${\sigma }_{{{\ell }}_{\mu }}({k}_{\mu })={\rm{arg}}{\rm{\Gamma }}({{\ell }}_{\mu }+1-{iZ}/{k}_{\mu })$ is the Coulomb phase shift for an electron with kinetic momentum ${k}_{\mu }=\sqrt{2{\varepsilon }_{\mu }}$ and angular momentum ${{\ell }}_{\mu }$ (${{\ell }}_{\alpha }$) escaping from the molecular ion with total charge Z. Asymptotically, the incoming wave ${{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}$ has outgoing spherical Coulomb waves only in channel μ, whereas there are incoming spherical Coulomb waves from all channels labelled with $\mu ^{\prime}$ = 1,..., N_c and ${{\mathcal{S}}}_{\mu {\mu }^{\prime }}^{{\rm{\Omega }}\dagger }={e}^{-2i{\delta }_{\mu }({k}_{\mu })}$ is the Hermitian conjugate of the scattering ${\mathcal{S}}$-matrix that contains the short range phase shift ${\delta }_{\mu }({k}_{\mu })$ (the procedure to obtain it is outlined below).

In contrast to atomic systems, for which evaluation of the electronic continuum lying just above the first ionization threshold is a single-channel scattering problem, in molecules it is always a multichannel scattering problem as a result of the non-central character of the molecular potential. For example, in the case of continuum states of ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ symmetry in the hydrogen molecule, there is an infinite number of open channels $\{\mu \}$ = [${{\rm{H}}}_{2}^{+}$($1s{\sigma }_{g}$), ${\ell }$ = 1, 3, 5, ..., $\infty$] in the energy region lying between the first ${{\rm{H}}}_{2}^{+}$($1s{\sigma }_{g}$) and the second ${{\rm{H}}}_{2}^{+}$($2p{\sigma }_{u}$) ionization thresholds (${E}_{1s{\sigma }_{g}}(R)\lt E\lt {E}_{2p{\sigma }_{u}}(R)$). For energies between the second and the third ionization thresholds, ${E}_{2p{\sigma }_{u}}(R)\lt E\lt {E}_{2p{\pi }_{u}}(R)$, new channels $\{\mu \}$ = [${{\rm{H}}}_{2}^{+}$($2p{\sigma }_{u}$), ${\ell }$ = 0, 2, 4, ..., $\infty$] open. The coupling between all these channels also prevents one from obtaining the corresponding states from a direct diagonalization of the electronic Hamiltonian ${\hat{{\mathcal{H}}}}_{{el}}$ in an ${{\mathcal{L}}}^{2}$ basis set, since the correct asymptotic boundary conditions of multichannel scattering states cannot be imposed in this way.

Our ${{\mathcal{L}}}^{2}$ multichannel solution first requires computation of the single-channel electronic uncoupled continuum states for each separate channel $\mu ^{\prime}$ in equation (27), the so-called uncoupled continuum states (UCSs). In other words, for each open ionic threshold α of H₂, a continuum state is evaluated for each value of angular momentum ${{\ell }}_{\mu }$. Although these states are unphysical, because of the lack of inter channel couplings, they can be renormalized easily to acquire the appropriate Dirac delta normalization of continuum states. In order to obtain scattering states ${\mathcal{P}}{{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}$ with appropriate incoming boundary conditions (28), the coupling between UCSs will then be introduced using the Lippmann–Schwinger formalism. More specifically, the UCSs are analogous to the static exchange wave function with continuum energy E

$$\begin{eqnarray}&&{\zeta }_{\mu E}^{0}={\rm{\Theta }}\left({\varphi }_{\mu }^{{\rm{\Omega }}}({{\bf{r}}}_{1},{\hat{{\bf{r}}}}_{2}){F}_{\mu E}^{{\rm{\Omega }}-}({r}_{2})\right).\end{eqnarray} \tag{ 30 }$$

For practical reasons, for each channel μ we first obtain the one-electron eigenstates of the ${{\ell }}_{\mu }$-projected Hamiltonian ${\hat{O}}_{{\hat{L}}^{2}}{\hat{h}}_{{el}}{\hat{O}}_{{\hat{L}}^{2}}$ (where ${\hat{O}}_{{\hat{L}}^{2}}$ is an angular momentum channel projection operator) instead of those from equation (16), thus reducing the multichannel problem to separately solving μ eigenvalue problems associated to the one-electron Hamiltonian diagonal blocks ${{\bf{h}}}_{i{{\ell }}_{\mu },i^{\prime} {{\ell }}_{\mu }}$ in equation (19). This allows us to obtain a basis set of ${{\ell }}_{\mu }$-projected molecular orbitals ${\left\{{\varphi }_{j}^{\lambda ,m,{{\ell }}_{\mu }}(r)\right\}}_{j=1}^{{N}_{B}}$ to be used as a basis to build two-electron configurations for a given single channel μ, i.e., ${\left\{{\rm{\Theta }}({\varphi }_{\mu }^{{\rm{\Omega }}}({{\bf{r}}}_{1},{\hat{{\bf{r}}}}_{2}){\varphi }_{j}({r}_{2}))\right\}}_{j=1}^{{N}_{B}}$(for simplicity in the notation, we have dropped the superscripts $\lambda ,m,{{\ell }}_{\mu }$ in ${\varphi }_{j}({r}_{2})$). The diagonalization of the two-electron Hamiltonian ${\hat{{\mathcal{H}}}}_{{el}}$ in terms of these static exchange configurations yields a discretized energy spectrum ${\{{E}_{\mu n}\}}_{n=1}^{{N}_{B}}$ and a set of discretized states (normalized to unity) in the form

$$\begin{eqnarray}&&{\tilde{\zeta }}_{\mu n}^{0}=\displaystyle \sum _{j=1}^{{N}_{B}}{C}_{j}^{\mu n}\left\{{\rm{\Theta }}({\varphi }_{\mu }^{{\rm{\Omega }}}({{\bf{r}}}_{1},{\hat{{\bf{r}}}}_{2}){\varphi }_{j}({r}_{2}))\right\}.\end{eqnarray} \tag{ 31 }$$

Since the channel label μ is fixed, the radial continuum wave function (discretized and normalized to unity) for electron 2 has the simple form ${\tilde{\xi }}_{\mu n}({r}_{2})$ = ${\displaystyle \sum }_{j}{C}_{j}^{\mu n}{\varphi }_{j}({r}_{2})$. Then the properly normalized discretized two-electron UCS with energy ${E}_{\mu n}$ state is obtained as

$$\begin{eqnarray}&&{\zeta }_{\mu {E}_{\mu n}}^{0}={\rho }_{\mu }^{1/2}({E}_{\mu n})\cdot {\rm{\Theta }}({\varphi }_{\mu }^{{\rm{\Omega }}}({{\bf{r}}}_{1},{\hat{{\bf{r}}}}_{2}){\tilde{\xi }}_{\mu n}({r}_{2})),\end{eqnarray} \tag{ 32 }$$

where the density of states ${\rho }_{\mu }({E}_{\mu n})$ = ${| \partial {E}_{\mu n^{\prime} }/\partial n^{\prime} | }_{n^{\prime} =n}^{-1}$ can be evaluated in most cases with a simple two-point interpolation formula ${\rho }_{\mu }({E}_{\mu n})$ = $2/({E}_{\mu (n+1)}-{E}_{\mu (n-1)})$. ${\zeta }_{\mu {E}_{\mu n}}^{0}$ is the ${{\mathcal{L}}}^{2}$ discretized approximation to the exact UCS ${\zeta }_{\mu E}^{0}$ in equation (30) within a box of size $r\in [0,{r}_{\mathrm{max}}]$. Accordingly, the radial part of the continuum wave function ${\xi }_{\mu {E}_{\mu n}}(r)$ = ${\rho }_{\mu }^{1/2}({E}_{\mu n})\cdot {\tilde{\xi }}_{\mu ,n}(r)$ is an approximation to the function ${F}_{\mu E}^{{\rm{\Omega }}-}({r}_{2})$. It should thus satisfy the boundary conditions given in equation (28). However, when calculations are performed in real space (as in most previously reported applications), the asymptotic form of the UCS continuum orbitals ${\xi }_{\mu {E}_{\mu n}}$ is instead

$$\begin{eqnarray}&&\underset{r\to \infty }{\mathrm{lim}}{\xi }_{\mu {E}_{\mu n}}(r)=\sqrt{\displaystyle \frac{2}{\pi {k}_{\mu n}}}\displaystyle \frac{1}{r}\mathrm{sin}\left[{\theta }_{\mu n}({k}_{\mu n})+{\delta }_{\mu n}({k}_{\mu n})\right],\end{eqnarray} \tag{ 33 }$$

where we are using ${E}_{\mu n}$-${E}_{\alpha }(R)$ = ${k}_{\mu n}^{2}/2$ for the excess kinetic energy of the ionized electron above the threshold α. In order to recover the appropriate complex asymptotic behavior of (28), one has to extract, from each UCS, the scattering phase shifts ${\delta }_{\mu n}$ associated with the short-range potential induced by the electrons that remain in the molecular ion and, with them available, the real continuum wave functions must be multiplied by the factor ${e}^{-i{\delta }_{\mu n}}$. If one is interested in obtaining observables such as angular distributions of the ejected electron, in which coherences between different partial waves manifest, the introduction of these phase shifts is crucial.

As we mentioned above, in the multichannel problem, the scattering states with the correct boundary conditions and energy E are related to the UCSs through the Lippmann–Schwinger equation:

$$\begin{eqnarray}&&{\mathcal{P}}{{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}={\zeta }_{\mu E}^{0}+{G}_{P}^{-}(E)V{\zeta }_{\mu E}^{0},\end{eqnarray} \tag{ 34 }$$

where the Green's function ${G}_{P}^{-}(E)$ must be calculated at the fixed energy E. In addition, in order to include the nuclear motion later, it is of utmost importance to have the same discretization $\{{E}_{\mu n}\}$ of the electron energy spectrum at all internuclear distances, i.e., ${E}_{\mu n}={E}_{\alpha }(R)+{\varepsilon }_{\mu n}$ for all R. In general, diagonalization of the UCS eigenvalue problem for a given μ channel does not fulfill the latter condition, unless by accident. To avoid this problem, we predefine a grid of electron kinetic energies, ${\{{\varepsilon }_{\mu i}\}}_{i=1}^{{N}_{\varepsilon }}$, and interpolate the UCS discretized continuum energies (as well as the shifts ${\delta }_{\mu n}$ and dipole transition moments between states) for each μ channel and each R. The latter implies that the box used to represent the system cannot be exactly the same for all μ and R. Indeed, every single-electron radial continuum function obtained by diagonalization inside the electronic box satisfies the boundary condition $r{\xi }_{\mu {E}_{\mu n}}(r={r}_{\mathrm{max}})$=0, whereas the analytical formula (33) with the interpolated phase shift ${\delta }_{\mu i}$ does not. Instead, the interpolated state shows its last node at a given radial distance r₀ very close to the edge of the box ${r}_{\mathrm{max}}$. To match one of the eigenvalues to a preselected value ${E}_{\mu i}$, the box length should be replaced by r₀. This would imply a different box size ${r}_{\mathrm{max}}$ and consequently a different basis set for each channel μ and each R. However, this is avoided with the introduction of a step potential ${V}_{0}{\rm{\Theta }}(r-{r}_{0})$ in the one-electron Hamiltonian ${\hat{h}}_{{el}}$ (${\rm{\Theta }}(x)$ is the Heaviside function and V₀ is a large real number), which guarantees that in the diagonalization procedure at least one of the eigenvalues of the new set of UCSs ${\zeta }_{\mu {E}_{\mu n}}^{0}$ matches the selected ${E}_{\mu i}$ value, but keeping the same box length and the same basis set for all μ and R.

A technical requirement to implement this interpolation method is that the density of states resulting from the direct diagonalization $\{{E}_{\mu n}\}$ must be similar or higher than that of the artificial energy grid $\{{E}_{\mu i}\}$. This condition prevents the same UCS ${\zeta }_{\mu {E}_{\mu n}}^{0}$ from being adjusted to several energy values ${E}_{\mu i}$, which may occasionally induce large displacements of the r₀ value from the box edge ${r}_{\mathrm{max}}$, with the subsequent loss of orthogonality among the continuum states. A common procedure in previous works is to use an ${N}_{\varepsilon }$-point energy grid ${\left\{{E}_{\mu i}={k}_{\mu i}^{2}/2+{E}_{\alpha }(R)\right\}}_{i=1}^{{N}_{\varepsilon }}$ with a linear spacing in momentum ${k}_{i}=i{\rm{\Delta }}k$, e.g., with ${\rm{\Delta }}k$ = 0.1 a.u. The latter choice (common to all angular momentum ${{\ell }}_{\mu }$ channel values and all ionization thresholds ${E}_{\alpha }(R)$ of the molecular ion) leads to a set of UCSs with an appropriate energy spacing for most applications using atto- and few-femtosecond XUV laser pulses with central frequencies ranging from ${\omega }_{c}$ = 16 eV up to 50 eV.

Finally, once the UCSs ${\tilde{\zeta }}_{\mu n}^{0}$ for all open channels $\{\mu \}$ are obtained, the coupled continuum wave function $P{{\rm{\Psi }}}_{\mu E}^{{\rm{\Omega }}-}$ for a given open channel μ is obtained by expressing the Lippmann–Schwinger equation (34) in terms of the UCS states, which takes the form

$$\begin{eqnarray}P{{\rm{\Psi }}}_{\mu {E}_{\mu n}}^{{\rm{\Omega }}-} & = & {\rho }_{\mu }^{1/2}({E}_{\mu n})\left({\tilde{\zeta }}_{\mu n}^{0}+\displaystyle \sum _{{\mu }^{\prime }n^{\prime} }\displaystyle \sum _{{\mu }^{\prime\prime }n^{\prime \prime} }\langle {\tilde{\zeta }}_{{\mu }^{\prime }n^{\prime} }^{0}| \right.{\hat{G}}_{P}^{-}({E}_{\mu n})\\ & & \times \left.| {\tilde{\zeta }}_{{\mu }^{\prime\prime }n^{\prime \prime} }^{0}\rangle \langle {\tilde{\zeta }}_{{\mu }^{\prime\prime }n^{\prime \prime} }^{0}| \hat{V}| {\tilde{\zeta }}_{\mu n}^{0}\rangle {\tilde{\zeta }}_{{\mu }^{\prime }n^{\prime} }^{0}\right),\end{eqnarray} \tag{ 35 }$$

where ${\hat{G}}_{P}^{-}(E)$ is the coupled Green's operator ${\hat{G}}_{P}^{-}(E)$ = ${\hat{G}}_{0}^{-}(E)$+${\hat{G}}_{0}^{-}(E)\hat{V}{\hat{G}}_{P}^{-}(E)$, which is written in terms of the Green's operator ${\hat{G}}_{0}^{-}(E)$ for the uncoupled channels, and $\hat{V}$ is the operator for the interchannel coupling (the non-block diagonal terms in the Hamiltonian ${\hat{{\mathcal{H}}}}_{{el}}$ connecting different channels $\mu \ne \mu ^{\prime}$). Detailed expressions for the Green's operator and the ensuing algebraic computational method can be found in [112, 114, 115]. Note that ${\rho }_{\mu }^{1/2}({E}_{\mu n})$ is introduced as defined in equation (32) to ensure the proper Dirac delta normalization of the continuum states.

As an illustration of the configurations typically used to describe one-photon processes, the UCS required to build the two lowest ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ continua associated to the ionic thresholds $1s{\sigma }_{g}\equiv 1{\sigma }_{g}\;\mathrm{and}\;2p{\sigma }_{u}\equiv 1{\sigma }_{u}$ of ${{\rm{H}}}_{2}$ are

• ${}^{1}{{\rm{\Sigma }}}_{u}^{+}[{H}_{2}^{+}(1s{\sigma }_{g}),\epsilon {\ell }{\text{}}]1{\sigma }_{g}n{\sigma }_{u}({\rm{n}}=1,75)$ up to four partial waves ${\ell }$ = 1, 3, 5, and 7.
• ${}^{1}{{\rm{\Sigma }}}_{u}^{+}[{H}_{2}^{+}(2p{\sigma }_{u}),\epsilon {\ell }{\text{}}]1{\sigma }_{u}n{\sigma }_{g}(n=1,75)$ up to four partial waves ${\ell }$ = 0, 2, 4, and 6.

In addition to the above, one has to calculate the dipole matrix elements between bound and continuum states for each $R$. The couplings between bound and continuum states ${\mathcal{P}}{{\rm{\Psi }}}_{\mu {E}_{\mu n}}^{{\rm{\Omega }}-}$ must be computed for each discretized continuum energy, i.e. for each of the ${N}_{\varepsilon }$ values in the set ${\{{E}_{\mu i}\}}_{i=1}^{{N}_{\varepsilon }}$. Moreover, to investigate above threshold ionization problems, where there are photon absorptions within the electronic continuum, it is also required to compute the continuum–continuum matrix elements $\langle {\mathcal{P}}{{\rm{\Psi }}}_{\mu {E}_{\mu i}}^{{\rm{\Omega }}-}| \hat{\mu }| {\mathcal{P}}{{\rm{\Psi }}}_{\mu {E}_{{\mu }^{\prime }j}}^{{\rm{\Omega }}-}\rangle$ in the double set of grid energies ${\{{E}_{\mu i},{E}_{{\mu }^{\prime }j}\}}_{i,j}^{{N}_{\varepsilon },{N}_{\varepsilon }^{\prime }}$, taking special care in the neighborhood of the points that satisfy ${E}_{\mu i}={E}_{{\mu }^{\prime }j}$, where a sufficiently large number of grid points must be used to account for the rapid variation of the corresponding dipole matrix elements with energy.

One of the most accurate and general methods to obtain eigenvalues and eigenstates in many-electron molecular systems is the configuration interaction (CI) procedure. We directly use antisymmetrized products of one-electron diatomic molecular (OEDM) wave functions, i.e. eigenstates of the parent molecular ion ${{\rm{H}}}_{2}^{+}$. Each antisymmetrized configuration in the total wavefunction ${{\rm{\Psi }}}^{{\rm{\Lambda }}}={\displaystyle \sum }_{{m}_{1},{m}_{2}/| {m}_{1}+{m}_{2}| ={\rm{\Lambda }}}$ ${\displaystyle \sum }_{{n}_{1}{n}_{2}}{C}_{{n}_{1}{n}_{2}}^{{m}_{1}{m}_{2}}{{\rm{\Psi }}}_{{n}_{1},{n}_{2}}^{{\rm{\Lambda }};{m}_{1},{m}_{2}}$ can be written in the form

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{n}_{1},{n}_{2}}^{{\rm{\Lambda }};{m}_{1},{m}_{2}}({{\bf{x}}}_{1},{{\bf{x}}}_{2})=\frac{1}{\sqrt{2}}\left\{{\varphi }_{{n}_{1}}^{{m}_{1}}({{\bf{r}}}_{1}){\varphi }_{{n}_{2}}^{{m}_{2}}({{\bf{r}}}_{2})\pm {\varphi }_{{n}_{2}}^{{m}_{2}}({{\bf{r}}}_{1}){\varphi }_{{n}_{1}}^{{m}_{1}}({{\bf{r}}}_{2})\right\}\\ &&\quad \times \frac{1}{\sqrt{2}}\left\{\alpha (1)\beta (2)\mp \beta (1)\alpha (2)\right\},\end{eqnarray} \tag{ 36 }$$

with ${\bf{x}}=\{{\bf{r}},\sigma \}$ for spatial ${\bf{r}}$ and spin $\sigma$ coordinates and ${\rm{\Lambda }}=| {m}_{1}+{m}_{2}|$ corresponds to the eigenvalue of the $z-{\text{}}{\rm{component}}\;{{\bf{L}}}_{z}$ of the total electronic angular momentum ${\bf{L}}={{\bf{l}}}_{1}+{{\bf{l}}}_{2}$ (with the known labels ${\rm{\Sigma }}({\rm{\Lambda }}=0)$, ${\rm{\Pi }}({\rm{\Lambda }}=1)$, ${\rm{\Delta }}({\rm{\Lambda }}=2)$,...). In the CI solution of a general two-electron problem, the construction of the total Hamiltonian only requires the computation of the electron–electron Coulomb repulsion. For the ${{\rm{H}}}_{2}$ molecule, we can reduce the problem, taking into account that the molecular configurations must also satisfy the point group symmetry requirements. For instance, states with ${{\rm{\Sigma }}}_{g}$ total symmetry require configurations of the type $[{\sigma }_{g}(1),{\sigma }_{g}(2)],[{\sigma }_{u}(1),{\sigma }_{u}(2)]$, $[{\pi }_{g}^{+1}(1),{\pi }_{g}^{-1}(2)]$, $[{\delta }_{g}^{+2}(1),{\delta }_{g}^{-2}(2)]$, etc. Similarly, to calculate states with ${{\rm{\Pi }}}_{u}^{+1}$ total symmetry the required configurations are $[{\sigma }_{g}(1),{\pi }_{u}^{+1}(2)],[{\sigma }_{u}(1),{\pi }_{g}^{+1}(2)]$, $[{\pi }_{u}^{-1}(1),{\delta }_{g}^{+2}(2)]$, $[{\pi }_{g}^{-1}(1),{\delta }_{u}^{+2}(2)]$, etc. Note that the computation of the matrix elements involving configurations of the type $[{\lambda }_{g/u}^{+{m}_{1}/-{m}_{1}},{\lambda }_{u/g}^{+{m}_{2}/-{m}_{2}}]{\text{}}{\rm{with}}\;\lambda \;\gt$ 0, four terms appear in the spatial part of equation (36) instead of two. Since we use a one-center expansion in partial waves for the OEDM representation, equation (17), the calculation of the two-electron integrals (direct plus exchange) $({\varphi }_{{n}_{1}}^{{\lambda }_{1},{m}_{1}}{\varphi }_{{n}_{2}}^{{\lambda }_{2},{m}_{2}}| 1/{r}_{12}| {\varphi }_{{n}_{3}}^{{\lambda }_{3},{m}_{3}}{\varphi }_{{n}_{4}}^{{\lambda }_{4},{m}_{4}})$ is analogous to that of the atomic case, by also using the expansion $1/{r}_{12}={\displaystyle \sum }_{k}({r}_{\lt }^{k}/{r}_{\gt }^{k+1}){P}_{k}(\mathrm{cos}{\theta }_{12})$. The radial part of the two electron integrals $\int {{\rm{d}}{r}}_{1}{{\rm{d}}{r}}_{2}{U}_{{n}_{1}{{\ell }}_{1}}({r}_{1}){U}_{{n}_{2}{{\ell }}_{2}}({r}_{2})({r}_{\lt }^{k}/{r}_{\gt }^{k+1}){U}_{{n}_{3}{{\ell }}_{3}}({r}_{1})$ ${U}_{{n}_{4}{{\ell }}_{4}}({r}_{2})$ is solved by first computing the function ${Y}^{k}({n}_{1}{{\ell }}_{1},{n}_{3}{{\ell }}_{3};{r}_{1})$ using the method of the auxiliary function ${Z}^{k}$ as in the standard method for atoms described in [116]. All angular integrations are analytical. This term involves a quintuple sum over ${{\ell }}_{1},{{\ell }}_{2},{{\ell }}_{3},{{\ell }}_{4}$ and k in addition to the typical number of two-electron integrals [${\mathcal{O}}({K}^{4}/8)$] (with $K$ the number of basis functions in which the molecular orbitals are expanded), which makes the evaluation of integrals very demanding.

The calculation of bound and resonant doubly excited states in our ${{\mathcal{L}}}^{2}$ procedure for ${{\rm{H}}}_{2}$ follows an analogous CI procedure, but the ${\mathcal{Q}}$ projection is performed by truncating the diagonalization to eliminate selected configurations. The doubly excited states, ${\psi }_{r}\equiv {\mathcal{Q}}{\psi }_{r}$ are indeed solutions of the Feshbach ${\mathcal{Q}}$-projected electronic Hamiltonian

$$\begin{eqnarray}&&\;[{\mathcal{Q}}{{\mathcal{H}}}_{{el}}{\mathcal{Q}}-{E}_{r}]{\psi }_{r}({{\bf{r}}}_{1},{{\bf{r}}}_{2};R)=0.\end{eqnarray} \tag{ 37 }$$

The completeness of the ${\mathcal{Q}}$ subspace is guaranteed by the projector ${\mathcal{Q}}={\displaystyle \sum }_{r^{\prime} }| {\psi }_{r^{\prime} }\rangle \langle {\psi }_{r^{\prime} }|$, built from all ${\mathcal{Q}}$ eigenvectors. In practice, to construct a single ${\mathcal{Q}}$ resonant space can be rather cumbersome. Instead, it is easier to partition ${\mathcal{Q}}$ and calculate separately each Rydberg series of resonant states ${\{{\psi }_{r}\}}_{{{\mathcal{Q}}}_{\alpha }}$ lying energetically above the ionization threshold ${E}_{\alpha }(R)$ and below the next one ${E}_{\alpha +1}(R)$. For a two-electron system, the Feshbach ${\mathcal{P}}$ projection operator for threshold $\alpha ,{{\mathcal{P}}}_{\alpha }(1,2)\;=$ $1-{{\mathcal{Q}}}_{\alpha }(1,2)$ can also be written in terms of one-particle projection operators, i.e., ${{\mathcal{P}}}_{\alpha }(1,2)\;=$ ${{\mathcal{P}}}_{\alpha }(1)+{{\mathcal{P}}}_{\alpha }(2)-{{\mathcal{P}}}_{\alpha }(1){{\mathcal{P}}}_{\alpha }(2)$, with ${{\mathcal{P}}}_{\alpha }(i)=| {\varphi }_{\alpha }(i)\rangle \langle {\varphi }_{\alpha }(i)|$ and where ${\varphi }_{\alpha }(i)$ is the molecular orbital associated to the ionic channel $\alpha$. Then molecular resonant states are built up in a hierarchy of resonance subspaces, namely, the set of ${{\mathcal{Q}}}_{\alpha =1}$ resonances appear by using the projector ${{\mathcal{P}}}_{\alpha =1}(i)=| 1s{\sigma }_{g}(i)\rangle \langle 1s{\sigma }_{g}(i)|$, then ${{\mathcal{Q}}}_{\alpha =2}$ resonance series with the projector ${{\mathcal{P}}}_{\alpha =1}(i)+{{\mathcal{P}}}_{\alpha =2}(i)=| 1s{\sigma }_{g}(i)\rangle \langle 1s{\sigma }_{g}(i)|$ $+\;| 2p{\sigma }_{u}(i)\rangle \langle 2p{\sigma }_{u}(i)|$ and, in general, ${{\mathcal{Q}}}_{\alpha =n}$ with the operator ${\displaystyle \sum }_{n}{{\mathcal{P}}}_{\alpha =n}(i)$. For instance, the lowest two series ${{\mathcal{Q}}}_{1}\;\mathrm{and}\;{{\mathcal{Q}}}_{2}{\text{}}{\;{\rm{of}}\;}^{1}{{\rm{\Sigma }}}_{u}^{+}$ resonances have been calculated using the following configurations

• ${{\mathcal{Q}}}_{1}^{1}{{\rm{\Sigma }}}_{u}^{+}$: (${{\mathcal{Q}}}_{1}$ projection removes all configurations that contain the $1{\sigma }_{g}$ orbital) $1{\sigma }_{u}n{\sigma }_{g}(n=2,70)$, $1{\pi }_{u}n{\pi }_{g}(n=1,70)$, $2{\sigma }_{g}n{\sigma }_{u}(n=2,35)$, $2{\sigma }_{u}n{\sigma }_{g}(n=2,18)$, $3{\sigma }_{g}n{\sigma }_{u}$ (n = 1, 18), $1{\pi }_{g}n{\pi }_{u}(n=1,10)$, $1{\delta }_{g}n{\delta }_{u}(n=1,10)$ and $2{\pi }_{u}n{\pi }_{g}(n=1,10)$.
• ${{\mathcal{Q}}}_{2}^{1}{{\rm{\Sigma }}}_{u}^{+}$: (${{\mathcal{Q}}}_{2}$ projection removes all configurations that contain the $1{\sigma }_{g}$ and/or the $1{\sigma }_{u}$ orbitals) $1{\pi }_{u}n{\pi }_{g}(n=1,70)$, $2{\sigma }_{g}n{\sigma }_{u}(n=2,70)$, $2{\sigma }_{u}n{\sigma }_{g}(n=3,18)$, $3{\sigma }_{g}n{\sigma }_{u}$ (n = 3, 18), $1{\pi }_{g}n{\pi }_{u}(n=2,10)$, $1{\delta }_{g}n{\delta }_{u}(n=1,10)$ and $2{\pi }_{u}n{\pi }_{g}(n=2,10)$.

It turns out that, by construction, the ${{\mathcal{Q}}}_{\alpha }$ space is strictly orthogonal to the ${{\mathcal{P}}}_{\alpha ^{\prime} }{\text{}}\;{\rm{when}}\;\alpha \geqslant \alpha ^{\prime}$ but they are not if $\alpha \lt \alpha ^{\prime}$, unless they correspond to subspaces of different molecular symmetry ${\rm{\Omega }}\;\mathrm{and}\;{\rm{\Omega }}^{\prime}$, for which orthogonality is automatically satisfied. As an illustration, the doubly excited states ${{\mathcal{Q}}}_{2}^{1}{{\rm{\Sigma }}}_{u}^{+}$ are orthogonal to the continua ${{\mathcal{P}}}_{1}{\ }^{1}{{\rm{\Sigma }}}_{u}^{+}(1s{\sigma }_{g},\epsilon {{\ell }}_{u})$ and ${{\mathcal{P}}}_{2}{\ }^{1}{{\rm{\Sigma }}}_{u}^{+}(2p{\sigma }_{u},\epsilon {{\ell }}_{g})$, but the doubly excited states ${{\mathcal{Q}}}_{1}^{\;1}{{\rm{\Sigma }}}_{u}^{+}$ are not strictly orthogonal to the ${{\mathcal{P}}}_{2}$ continuum space. Therefore, rigorously, the corresponding overlaps should be included in the solution of the system of coupled differential equations (13), although in practice such overlaps are usually small and can often be neglected. According to the formal Feshbach theory, the ${\mathcal{Q}}$ subspace contains only the resonance contributions and the complementary space ${\mathcal{P}}$ all the rest, which means that molecular bound states pertain to the ${\mathcal{P}}$ space. However, the latter states are of poorer quality than those resulting from a direct diagonalization of the unprojected electronic Hamiltonian, since configurations representing double excitations can be important to correctly describe the ground state as well as the Rydberg series of singly excited states. Of course, since bound states constructed in this way are not exactly orthogonal to either the ${{\mathcal{Q}}}_{\alpha }$ or the ${{\mathcal{P}}}_{\alpha }$ subspaces, the corresponding overlaps should also be included (although, as for the states associated to different ${{\mathcal{Q}}}_{\alpha }$ subspaces, the overlaps are small and may be neglected). In homonuclear diatomic molecules subject to one-photon processes, the selection rule g $\leftrightarrow$ u for photoionization from any bound state formally eliminates the problem. However, in multiphoton absorption processes and laser XUV-pump-IR-probe schemes the consequences of neglecting the overlaps must be tested in each particular case.

The expansion of the time-dependent wave packet, equation (12), in a truncated set of vibronic states, $\{{\psi }_{n}({\bf{r}};R)\cdot {\chi }_{{v}_{n}}(R)\}$, requires the computation of a complete set of nuclear wave functions ${\chi }_{{v}_{n}}$ for each (bound, resonant and continuum) electronic state. This comprises both bound and dissociating vibrational states associated to the electronic potential energy curves (PEC) of the neutral and singly charged diatomic molecule. The vibrational states are thus obtained by solving, for each potential energy surface ${E}_{n}(R)$ and rotational number $J$, the single channel equation

$$\begin{eqnarray}&&\left[\frac{{d}^{2}}{{{dR}}^{2}}-\frac{J(J+1)}{{R}^{2}}-2\mu ({E}_{n}(R)-{W}_{{{nv}}_{n}})\right]{\chi }_{{v}_{n}}(R)=0,\end{eqnarray} \tag{ 38 }$$

where $\mu$ is the nuclear reduced mass, ${W}_{{{nv}}_{n}}$ is the total vibronic energy and ${E}_{n}(R)$ is the $n$ th electronic PEC. According to our previous notation in equation (12), for vibrational states supported by electronic bound states we use the notation $n=b$, for vibrational states associated to doubly excited states, $n=r$, and for vibrational states associated to each ionic state, $n=\alpha$.

In order to solve the nuclear equation, the wave functions ${\chi }_{{v}_{n}}(R)$ are expanded in terms of ${N}_{R}$ B-splines of order $k$, ${\chi }_{{v}_{n}}(R)={\displaystyle \sum }_{i=1}^{{N}_{R}}{c}_{i}^{{v}_{n}}{B}_{i}^{k}(R)$, and equation (38) is transformed into a matrix eigenvalue problem, which is readily diagonalized. Typical numbers used in previous calculations are ${N}_{R}$ = 240 B-splines of order $k$ = 8 inside a box of ${R}_{\mathrm{max}}$ = 14 a.u. Eigenvalues lying above the dissociation limit of each PEC correspond to discretized continuum radial states ${\tilde{\chi }}_{{v}_{n}}(R)$, where the tilde indicates that they are normalized to the Kronecker delta, because they come from a diagonalization using an ${{\mathcal{L}}}^{2}$ basis. To recover the proper Dirac delta normalization in this effective one-channel nuclear problem, the discretized continuum states are then renormalized through the density of states, i.e., ${\chi }_{{v}_{n}}={\rho }_{{v}_{n}}^{1/2}{\tilde{\chi }}_{{v}_{n}}\;{\rm{where}}\;{\rho }_{{v}_{n}}$ is the density of continuum states given by ${\rho }_{{v}_{n}}=2/({W}_{n({v}_{n}+1)}-{W}_{n({v}_{n}-1)})$. These continuum radial wave functions must also satisfy asymptotic incoming boundary conditions (photodissociation must be understood as a half collision, without entrance channels) in the form

$$\begin{eqnarray}&&\underset{R\to \infty }{\mathrm{lim}}{\chi }_{{v}_{n}}(R)=\sqrt{\frac{\mu }{2\pi {K}_{n}}}\frac{1}{{iR}}\left({e}^{i{\theta }_{n}}-{{\mathcal{S}}}_{n}^{\dagger }{e}^{-i{\theta }_{n}}\right),\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}&&{\theta }_{n}={K}_{n}R-\frac{J\pi }{2}.\end{eqnarray} \tag{ 40 }$$

${K}_{n}$ is the nuclear momentum and the ${\mathcal{S}}$-matrix contains the short-range nuclear phase shift ${{\rm{\Delta }}}_{n}$, ${{\mathcal{S}}}_{n}=\mathrm{exp}(-2i{{\rm{\Delta }}}_{n}({K}_{n}))$, which can be explicitly calculated by adjusting the real nuclear radial wave function resulting from the diagonalization to the known analytical asymptotic form

$$\begin{eqnarray}&&\underset{R\to \infty }{\mathrm{lim}}{\chi }_{{v}_{n}}(R)=\sqrt{\frac{2\mu }{\pi {K}_{n}}}\frac{1}{R}\mathrm{sin}\left[{\theta }_{n}({K}_{n})+{{\rm{\Delta }}}_{n}({K}_{n})\right]\end{eqnarray} \tag{ 41 }$$

in the outer part of the nuclear radial box. Due to the incoming boundary conditions implicit in the ${\mathcal{S}}$-matrix, our real vibrational wave functions must then be multiplied by the extra phase factor ${e}^{-i{{\rm{\Delta }}}_{n}}$ to satisfy equation (39).

3.6.1. Box size

The interaction time with the electromagnetic field largely determines the spatial dimension required to describe the electron and nuclear dynamics in numerical methods using finite size boxes. For the variety of open channels in ${{\rm{H}}}_{2}$ photoionization,

$$\begin{eqnarray*}&&\begin{array}{l}{\rm{direct}}\;{\rm{dissociative}}\\ \quad {\rm{ionization}}\end{array} & {{\rm{H}}}_{2}+{\rm{n}}{\rm{\hslash }}\omega \to {\rm{H}}{+H}^{+}{+e}^{-}\\ &&\begin{array}{l}{\rm{direct}}\;{\rm{non-dissociative}}\\ \quad \;{\rm{ionization}}\end{array} & {{\rm{H}}}_{2}+{\rm{n}}{\rm{\hslash }}\omega \to {{\rm{H}}}_{2}^{+}{+e}^{-}\\ {\rm{autoionization}} & {{\rm{H}}}_{2}+{\rm{n}}{\rm{\hslash }}\omega \to {{\rm{H}}}_{2}^{**}\to {\rm{H}}{+{\rm{H}}}^{+}{+e}^{-}\\ & {{\rm{H}}}_{2}+{\rm{n}}{\rm{\hslash }}\omega \to {{\rm{H}}}_{2}^{**}\to {{\rm{H}}}_{2}^{+}{+{\rm{e}}}^{-}\end{eqnarray*}$$

and even in dissociative excitation through the singly or doubly excited states (${{\rm{H}}}_{2}$+ n ${\rm{\hslash }}\omega \to {{\rm{H}}}_{2}^{**}\to$ H + H or ${{\rm{H}}}_{2}$+ $n\;{\rm{\hslash }}\omega \to {{\rm{H}}}_{2}^{**}\to {{\rm{H}}}^{+}$ +H⁻; ${{\rm{H}}}_{2}$+ n ${\rm{\hslash }}\omega \to {{\rm{H}}}_{2}^{*}\to$ H + H), we have electron and/or nuclear fragments gaining a given kinetic energy, which makes their associated wave packets to travel in their respective electron and nuclear boxes (${r}_{\mathrm{max}}\;\mathrm{and}\;{R}_{\mathrm{max}}$). Long time propagations may produce unphysical reflections of these wave packets at the edge of their boxes. The energy condition, the counterpart of this time limitation, is that the energy spacing between discretized states (both electronic and nuclear) must be smaller than the spectral width ${\rm{\Delta }}\omega$ of the laser pulse given by its Fourier transform, as discussed in section 2. Similarly, the energy spacing between non-resonant ${\mathcal{P}}$ electronic continuum states must be smaller than the autoionization widths of the relevant doubly excited states in the range of internuclear distances where autoionization occurs.

3.6.2. Gauge choice

3.6.3. Wave packet normalization

The TDSE must be solved by integrating the coupled equations (13) up to a time ${t}_{f}$ in which the stationary regime is reached. This time must be larger than the pulse duration $T$ and larger than the longest lifetime $1/{{\rm{\Gamma }}}^{\lt }$ of the populated doubly excited states (DES) (${t}_{f}\gt \{T,1/{{\rm{\Gamma }}}^{\lt }\}$). The physical information on the photoionization process is contained in the expansion coefficients ${C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}(t={t}_{f})$, which represent the transition amplitudes to the final continuum states and from which all relevant observables can be extracted, from total cross-sections to ionization probabilities, either integrated or differential in ion and/or electron kinetic energies as well as in ion and/or electron ejection angles.

3.7.1. Ionization probabilities

The probability of finding the molecule in a given asymptotic state at $t={t}_{f}$ is obtained by projecting the time-dependent wave function into that state. At $t={t}_{f}$, autoionizing states have either decayed completely into the ionization continua or have led to the dissociation of the molecule into two neutral fragments. In either case, the asymptotic state can be written, to a very good approximation, as $\{{\psi }_{n}({\bf{r}},R)\cdot {\chi }_{{v}_{n}}(R)\}$ (as defined in section 3.1). Because the time-dependent wave function (12) is already expanded in the basis of asymptotic states, the probability is given directly by the square of the expansion coefficient at $t={t}_{f}$. A practical way of choosing ${t}_{f}$ is to check that the projections of the molecular wave packet into the ${\mathcal{P}}\;\mathrm{and}\;{\mathcal{Q}}$ subspaces remain constant. This is equivalent to reaching the asymptotic limit, which in our Feshbach context is equivalent to ensuring that ${\mathcal{Q}}\hat{{\mathcal{H}}}{\mathcal{P}}+{\mathcal{P}}\hat{{\mathcal{H}}}{\mathcal{Q}}\to 0{\text{}}\;{\rm{ast}}\;\to \infty$. In this limit, the double differential ionization probability, differential in both vibrational (${E}_{{v}_{\alpha }}$) and electronic (${\varepsilon }_{\alpha }$) energies, is written, for each ${\ell }$-partial wave of the ionized electron, as

$$\begin{eqnarray}&&\frac{{d}^{2}{P}^{{{\ell }}_{\alpha }}({E}_{{v}_{\alpha }},{\varepsilon }_{\alpha },T)}{{{dE}}_{{v}_{\alpha }}d{\varepsilon }_{\alpha }}={| {C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}({t}_{f})| }^{2},\end{eqnarray} \tag{ 42 }$$

where ${C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}$ is the expansion coefficient that multiplies the asymtotic continuum state $\{{\psi }_{\alpha {{\ell }}_{\alpha }}^{{\varepsilon }_{\alpha }}\cdot {\chi }_{{v}_{\alpha }}\}$ with total energy ${W}_{\alpha {v}_{\alpha }}={E}_{{v}_{\alpha }}+{\varepsilon }_{\alpha }$.

The singly differential ionization probabilities are obtained by performing an additional integral/sum over the labels corresponding either to electronic or vibrational states in equation (42). For instance, the nuclear kinetic energy distributions are calculated from the expression:

$$\begin{eqnarray}&&\frac{{dP}({E}_{{v}_{\alpha }},T)}{{{dE}}_{{v}_{\alpha }}}=\displaystyle \sum _{{{\ell }}_{\alpha }}\int {\rm{d}}{\varepsilon }_{\alpha }{| {C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}({t}_{f})| }^{2},\end{eqnarray} \tag{ 43 }$$

where we have integrated over all electronic continuum energies ${\varepsilon }_{\alpha }$ and summed over all electron partial waves ${{\ell }}_{\alpha }$ included in the expansion (12). An analogous expression can be obtained for the electron kinetic energy distributions:

$$\begin{eqnarray}&&\frac{{dP}({\varepsilon }_{\alpha },T)}{d{\varepsilon }_{\alpha }}=\displaystyle \sum _{{{\ell }}_{\alpha }}{\displaystyle \sum \int }_{{v}_{\alpha }}{| {C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}({t}_{f})| }^{2},\end{eqnarray} \tag{ 44 }$$

where the symbol $\sum \int$ indicates a sum over the bound vibrational states and an energy integral over all scattering nuclear states lying above the dissociation limit.

Another interesting observable concerns the electron angular distributions in the photoionization process. The asymptotic wave function satisfying incoming wave boundary conditions that describes an ejected electron with energy ${\varepsilon }_{\alpha }$ in channel $\alpha$ and emitted in the direction ${\hat{{\bf{k}}}}_{e}\equiv d{{\rm{\Omega }}}_{e}$ can be written as a partial-wave decomposition in spherical harmonics in the form [118]

$$\begin{eqnarray}{{\rm{\Psi }}}_{\alpha }^{(-)}({\bf{r}},R,t) & = & \displaystyle \sum _{{{\ell }}_{\alpha }{m}_{\alpha }}{i}^{{{\ell }}_{\alpha }}{e}^{-i{\sigma }_{{{\ell }}_{\alpha }}({\varepsilon }_{\alpha })}{{\mathcal{Y}}}_{{{\ell }}_{\alpha }}^{{m}_{\alpha }*}({\hat{{\bf{k}}}}_{e})\\ & & \times {\psi }_{\alpha {{\ell }}_{\alpha }{m}_{\alpha }}^{{\varepsilon }_{\alpha }}({\bf{r}},R){\chi }_{{v}_{\alpha }}(R){e}^{-{{iW}}_{\alpha {v}_{\alpha }}t},\end{eqnarray} \tag{ 45 }$$

where ${\sigma }_{{{\ell }}_{\alpha }}$ is the Coulomb phase shift given after equation (29). Accordingly, to obtain the triply differential ionization probabilities, differential in the ejection solid angle $d{{\rm{\Omega }}}_{e}$ and electronic ${\varepsilon }_{\alpha }$ and nuclear energies ${E}_{{v}_{\alpha }}$, one has to project the time-dependent wave function (12) into the asymptotic expression (45) to yield

$$\begin{eqnarray}&&\frac{{d}^{3}P}{{{dE}}_{{v}_{\alpha }}d{\varepsilon }_{\alpha }d{{\rm{\Omega }}}_{e}}={\left|\displaystyle \sum _{\alpha ,{{\ell }}_{\alpha },{m}_{\alpha }}{i}^{-{{\ell }}_{\alpha }}{e}^{i{\sigma }_{{{\ell }}_{\alpha }}({\varepsilon }_{\alpha })}{{\mathcal{Y}}}_{{{\ell }}_{\alpha }}^{{m}_{\alpha }}({\hat{{\bf{k}}}}_{e}){C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}({t}_{f})\right|}^{2}.\end{eqnarray} \tag{ 46 }$$

This fully differential (in all angles and energies) probability provides the most detailed information on the intrinsic dynamics of the target [22, 86], information that can be hidden in the integrated probabilities. From an experimental point of view, it implies that all fragments must be detected in coincidence to retrieve their individual momenta. Microscopes making use of velocity correlation (VC) methods [119] or velocity map imaging (VMI) [120, 121] are able to provide the velocity vectors for a single event, e.g., for the ejection of an ion/electron pair upon photoionization of a single target molecule. Coincidence measurements are also possible through using reaction microscopes, such as those based on cold target recoil ion momentum xpectroscopy (COLTRIMS) [99, 122]. A peculiarity of COLTRIMS experiments is that one can distinguish between ions ejected, for example, in the left and right directions (or in the upward and downward directions, which is the usual convention we use for molecules oriented along the $Z$ axis). This breaks the inversion symmetry of homonuclear diatomic molecules. For example, in the case of ${{\rm{H}}}_{2}$ dissociative ionization, a COLTRIMS experiment allows one to know in which atomic center the residual electron is localized after dissociation. The two possible channels are H($n{\ell }$)+H⁺(right) and H⁺(left)+H($n{\ell }$). In order to obtain the correct ionization probabilities associated to these channels one has to project the time-dependent wave function in the combination of $g\;\mathrm{and}\;u$ asymptotic channels that lead to either of these channels. When only the $\alpha =1s{\sigma }_{g}\;\mathrm{and}\;\alpha =2p{\sigma }_{u}$ channels are accessible, the appropriate combination is

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{right,left}}}^{(-)}=\frac{1}{\sqrt{2}}[{{\rm{\Psi }}}_{1s{\sigma }_{g}}^{(-)}\pm {{\rm{\Psi }}}_{2p{\sigma }_{u}}^{(-)}],\end{eqnarray} \tag{ 47 }$$

which leads to H $(1{\rm{s}}){+{\rm{H}}}^{+}({\rm{right}})\;{\rm{and}}\;{{\rm{H}}}^{+}({\rm{left}})+{\rm{H}}(1{\rm{s}})$. This implies that the ionization amplitudes ${C}_{1s{\sigma }_{g}{{\ell }}_{1s{\sigma }_{g}}{v}_{1s{\sigma }_{g}}}^{{\varepsilon }_{1s{\sigma }_{g}}}(t)$ and ${C}_{2p{\sigma }_{u}{{\ell }}_{2p{\sigma }_{u}}{v}_{2p{\sigma }_{u}}}^{{\varepsilon }_{2p{\sigma }_{u}}}(t)$ must be also combined.

3.7.2. Cross-sections

For relatively long and not too intense pulses, it is sometimes useful to evaluate cross-sections rather than probabilities. The cross-section for single ionization due to an $N$-photon absorption process is given, in the framework of perturbation theory, by:

$$\begin{eqnarray}&&{\sigma }^{(N)}=\frac{{{\rm{\Gamma }}}^{(N)}}{{F}^{N}},\end{eqnarray} \tag{ 48 }$$

where $F$ stands for the photon flux per unit of area and time and ${{\rm{\Gamma }}}^{(N)}$ is the transition rate for a multiphoton process given by the Fermi's golden rule:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{(N)}=\frac{2\pi }{{\rm{\hslash }}}\left|\langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| \hat{\mu }\frac{1}{{E}_{i}+(N-1){\rm{\hslash }}\omega -\hat{{\mathcal{H}}}}\right.{\left.\hat{\mu }...\frac{1}{{E}_{i}+{\rm{\hslash }}\omega -\hat{{\mathcal{H}}}}\hat{\mu }| {{\rm{\Psi }}}_{{E}_{i}}\rangle \right|}^{2},\end{eqnarray} \tag{ 49 }$$

where ${{\rm{\Psi }}}_{{E}_{i}}$ is the initial (ground) vibronic state with total energy ${E}_{i}\equiv {W}_{{{gv}}_{g}}\;\mathrm{and}\;{{\rm{\Psi }}}_{{E}_{f}}^{-}$ corresponds to the (energy normalized) final vibronic state within the electronic continuum with energy ${E}_{f}\equiv {W}_{\alpha {v}_{\alpha }}$, and $\hat{\mu }$ is the dipole operator that can be written in the velocity form ${\hat{\mu }}_{V}=\displaystyle \frac{e}{{mc}}{\bf{p}}$ or in the length form ${\hat{\mu }}_{L}=e{\bf{r}}$.

The one-photon cross-section, which is given in units of length${}^{2}$, can be written in the velocity gauge as

$$\begin{eqnarray}&&\frac{d\sigma }{d{\rm{\Omega }}}=\frac{4{\pi }^{2}\alpha {\rm{\hslash }}}{{m}^{2}{\omega }_{{fi}}}{| \langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{p}}| {{\rm{\Psi }}}_{{E}_{i}}\rangle | }^{2},\end{eqnarray} \tag{ 50 }$$

and in the length gauge as

$$\begin{eqnarray}&&\frac{d\sigma }{d{\rm{\Omega }}}=4{\pi }^{2}\alpha {\rm{\hslash }}{\omega }_{{fi}}{| \langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{r}}| {{\rm{\Psi }}}_{{E}_{i}}\rangle | }^{2},\end{eqnarray} \tag{ 51 }$$

where ${\omega }_{{fi}}=({E}_{f}-{E}_{i})/{\rm{\hslash }},\alpha ={e}^{2}/{\rm{\hslash }}$c is the fine-structure constant and ${{\boldsymbol{\epsilon }}}_{p}$ is the laser field polarization vector. These expressions correspond to generalized cross-sections, i.e., they are independent of the field parameters. One should then be able to extract the dipole transition matrix elements to the continuum from the one-photon ionization amplitudes ${C}_{{E}_{f}}^{(N=1)}(t)[\equiv {C}_{\alpha {{\ell }}_{\alpha }{v}_{\alpha }}^{{\varepsilon }_{\alpha }}(t)]$ resulting from the solution of the time-dependent Schrödinger equation [94, 95] (see equation (12)). The superscript $(N=1)$ in the amplitudes specifies here a one-photon transition amplitude. For the sake of conciseness, here we only give the expressions for the velocity gauge. For a one-photon absorption process driven by a laser pulse of duration $T$ and carrier frequency ${\omega }_{c}$, represented with a vector potential ${\bf{A}}(t)={A}_{0}{{\mathfrak{F}}}_{{\omega }_{c}}(t){{\boldsymbol{\epsilon }}}_{p}$, one can relate the one-photon ionization amplitudes to the dipole transition matrix elements using first-order time-dependent perturbation theory through the expression

$$\begin{eqnarray}&&{P}_{{E}_{f}}({t}_{f})={| {C}_{{E}_{f}}^{(N=1)}({t}_{f})| }^{2}\\ &&={\left|\frac{-{{ieA}}_{0}}{{\rm{\hslash }}{mc}}\langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{p}}| {{\rm{\Psi }}}_{{E}_{i}}\rangle {K}^{(1)}({\omega }_{c},{\omega }_{{fi}},T)\right|}^{2},\end{eqnarray} \tag{ 52 }$$

in which the only time-dependent term is ${K}^{(1)}({\omega }_{c},{\omega }_{{fi}},T)$. The time ${t}_{f}(\gt T)$ corresponds to the final integration time previously introduced. The kernel ${K}^{(1)}$ [95] or shape function [94] corresponds to the ${\omega }_{{fi}}$-Fourier transform of the function ${{\mathfrak{F}}}_{{\omega }_{c}}(t)$ that enters in the definition of the vector potential ${\bf{A}}(t)$ defined in the interval $t\in [0,T]$, i.e.,

$$\begin{eqnarray}&&{K}^{(1)}({\omega }_{c},{\omega }_{{fi}},T)={\displaystyle \int }_{0}^{T}{e}^{i{\omega }_{{fi}}t}{{\mathfrak{F}}}_{{\omega }_{c}}(t){\rm{d}}{t}.\end{eqnarray} \tag{ 53 }$$

By substituting the expressions (52) and (53) in equation (50), the one-photon ionization cross-section can be expressed now in terms of the ionization probability as follows

$$\begin{eqnarray}&&\frac{d\sigma }{d{\rm{\Omega }}}=\frac{4{\pi }^{2}{e}^{2}{\rm{\hslash }}}{{\omega }_{{fi}}\alpha }\frac{{| {C}_{{E}_{f}}^{(N=1)}({t}_{f})| }^{2}}{{| {A}_{0}| }^{2}{| {K}^{(1)}({\omega }_{c},{\omega }_{{fi}},T)| }^{2}}.\end{eqnarray} \tag{ 54 }$$

In the limit of long pulse duration (cw radiation), the term ${| {K}^{(1)}({\omega }_{c},{\omega }_{{fi}},T)| }^{2}$ tends to $(3\pi T/16)\delta ({\omega }_{c}-{\omega }_{{fi}})$, in agreement with the long-time normalization coefficients used in [78]. Equivalently, from the expressions (48) and (49), the two-photon ionization cross-section in units of length${}^{4}\cdot$ time, can be written as

$$\begin{eqnarray}\frac{d\sigma }{d{\rm{\Omega }}} & = & 8{\pi }^{3}{\alpha }^{2}{{\rm{\hslash }}}^{3}{({\omega }_{{fi}}/2)}^{2}| \langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{r}}\\ & & \times {[{E}_{i}+{\rm{\hslash }}{\omega }_{c}-\hat{{\mathcal{H}}}]}^{-1}{{\boldsymbol{\epsilon }}}_{p}\cdot {{\bf{r}}| {{\rm{\Psi }}}_{{E}_{i}}\rangle | }^{2}\end{eqnarray} \tag{ 55 }$$

in the length gauge and

$$\begin{eqnarray}\frac{d\sigma }{d{\rm{\Omega }}} & = & \frac{8{\pi }^{3}{\alpha }^{2}{{\rm{\hslash }}}^{3}}{{m}^{4}{({\omega }_{{fi}}/2)}^{2}}| \langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{p}}\\ & & \times {[{E}_{i}+{\rm{\hslash }}{\omega }_{c}-\hat{{\mathcal{H}}}]}^{-1}{{\boldsymbol{\epsilon }}}_{p}\cdot {{\bf{p}}| {{\rm{\Psi }}}_{{E}_{i}}\rangle | }^{2}\end{eqnarray} \tag{ 56 }$$

in the velocity gauge. Again, to avoid unnecessary repetitions, we only consider the velocity gauge henceforth. The amplitudes extracted from a time-dependent calculation and the transition matrix elements given by the Fermi's golden rule are related through the following expression from second order time-dependent perturbation theory:

$$\begin{eqnarray}{P}_{{E}_{f}}({t}_{f}) & = & {| {C}_{{E}_{f}}^{(N=2)}({t}_{f})| }^{2}=\left|{\left(\frac{-i}{{\rm{\hslash }}}\right)}^{2}{\left(\frac{{{eA}}_{0}}{{mc}}\right)}^{2}\right.\\ & & \times \displaystyle \sum _{m}\langle {{\rm{\Psi }}}_{{E}_{f}}^{-}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{p}}| {{\rm{\Psi }}}_{{E}_{m}}\rangle \langle {{\rm{\Psi }}}_{{E}_{m}}| {{\boldsymbol{\epsilon }}}_{p}\cdot {\bf{p}}| {{\rm{\Psi }}}_{{E}_{i}}\rangle \\ & & \times {\left.{K}^{(2)}({\omega }_{c},{\omega }_{{fm}},{\omega }_{{mi}},T)\right|}^{2},\end{eqnarray} \tag{ 57 }$$

where

$$\begin{eqnarray}{K}^{(2)}({\omega }_{c},{\omega }_{{fm}},{\omega }_{{mi}},T) & = & {\displaystyle \int }_{0}^{T}{\rm{d}}{t}^{\prime} {e}^{i{\omega }_{{fm}}t^{\prime} }{{\mathfrak{F}}}_{{\omega }_{c}}(t^{\prime} )\\ & & \times {\displaystyle \int }_{0}^{{t}^{\prime }}{\rm{d}}{t}^{\prime \prime} {e}^{i{\omega }_{{mi}}t^{\prime \prime} }{{\mathfrak{F}}}_{{\omega }_{c}}(t^{\prime \prime} ){\rm{d}}{t}^{\prime \prime} .\end{eqnarray} \tag{ 58 }$$

In general, the fact that ${K}^{(2)}$ depends on the intermediate state energies, ${E}_{m}={\rm{\hslash }}{\omega }_{m}$, prevents one from factorizing this term. However, an expression equivalent to equation (54) can be obtained when no intermediate $| {{\rm{\Psi }}}_{{E}_{m}}\rangle$ states are resonantly populated and the pulse is long enough (the longer $T$ the shorter the energy bandwidth). When these conditions are fulfilled, the dependence of ${K}^{(2)}({\omega }_{c},{\omega }_{{fm}},{\omega }_{{mi}},T){\text{}}{\;{\rm{on}}\;{E}}_{m}$ can be safely ignored and an approximate analytical expression $K{^{\prime} }^{(2)}({\omega }_{c},{\omega }_{{fi}},T)$, similar to that given in [94], can be used. By using this approximation, the two-photon single ionization cross-section can be written as

$$\begin{eqnarray}&&\frac{d\sigma }{d{\rm{\Omega }}}=\frac{8{\pi }^{3}{{\rm{\hslash }}}^{3}{\alpha }^{2}}{{\bar{(}{\omega }_{{fi}}/2)}^{2}{m}^{4}}\frac{{| {C}_{{E}_{f}}^{(N=2)}({t}_{f})| }^{2}}{{\left(\frac{e| {A}_{0}| }{{mc}{\rm{\hslash }}}\right)}^{4}{| K{^{\prime} }^{(2)}({\omega }_{c},{\omega }_{{fi}},T)| }^{2}}.\end{eqnarray} \tag{ 59 }$$

In the limit of $T\to \infty$, the function ${| K{^{\prime} }^{(2)}({\omega }_{c},{\omega }_{{fi}},T)| }^{2}$ tends to $(70\pi T/{64}^{2})\delta ({\omega }_{{fi}}-2{\omega }_{c})$. The latter asymptotic expression was used previously with different pulse durations to build an effective cross-section for comparison purposes [69, 72]. However, strictly speaking, when considering the vibrational structure of a molecular target, the condition for the separability of equation (57) is never fulfilled in practice. For this reason, in the description of multiphoton processes induced by ultrashort laser pulses, one usually discusses the results in terms of the calculated probabilities, ${P}_{{E}_{f}}({t}_{f})={| {C}_{{E}_{f}}^{(N=2)}({t}_{f})| }^{2}$.

In 2009, UV-pump/IR-probe experiments were performed to track dissociative photoionization of ${{\rm{H}}}_{2}$ and D${}_{2}$ using femtosecond pulses, which allowed the investigation of the NWP dynamics associated to field-dressed ionic states [143]. One year later, the faster electron dynamics associated to excitation in the region of the D${}_{2}$ autoionizing states was explored in a similar two-color UV–IR scheme, but now using an isolated attosecond UV pulse as the pump [76]. While in the first problem the theoretical description can be achieved by simply considering the NWP moving in a one-dimensional PEC, an accurate description of the second one mandatorily requires treating both electronic and nuclear motion. As we discussed in detail in sections 4 and 5, autoionization from DES and nuclear motion occur in a similar timescale, and consequently electronic and nuclear motions must be coupled during the time propagation. However, the addition of an IR field considerably complicates the picture with respect to the physics shown in the two previous sections, often accompanied by a significant increase in the computational effort. Indeed, an electron excited to the continuum by the UV pulse can absorb IR photons (above threshold ionization), thus gaining both energy and angular momenta. The former implies the use of larger radial boxes and denser grids, and the latter a dramatic increase in the number of partial waves included in the angular expansion of the wave function. The introduction of a full set of dipole couplings between continuum states of different symmetries (for instance, ${}^{1}{{\rm{\Sigma }}}_{g}^{+}{\leftrightarrow }^{1}{{\rm{\Sigma }}}_{u}^{+}$ for linear polarization along the molecular axis), in order to represent these above-threshold IR transitions, comes with the cost of including double sets of electronic and vibrational continua for each ionization threshold of H${}_{2}^{+}(n{\ell }{\lambda }_{g/u})$. In spite of these additional complications, the huge improvement in computational capabilities in the last few years, e.g., by means of supercomputers, have made it possible to perform the first calculations considering molecular ionization by broadband UV pulses combined with IR fields.

In 2010, the combination of these theoretical efforts with a pioneering experiment that made use of an isolated UV attosecond pump pulse and an IR fs probe pulse demonstrated the possibility of controling the localization of the electronic charge within the molecule as it dissociates [76]. In this specific UV–IR scheme, a single attosecond pulse of 300--400 as centered at 30 eV, whose spectrum extends from 20 to 40 eV, was used to trigger a prompt ionization. Such a bandwidth covers a wide range of energies, thus leading to population of the first two ionization continua and the first series of ${{\mathcal{Q}}}_{1}$ DES (see the energy diagrams plotted in figure 14). In the absence of the IR radiation, the angle- and energy-differential ionization probabilities follow the patterns discussed in section 4. However, when the UV pulse is combined with a time-delayed intense IR pulse with the same linear polarization and a FWHM of 6 fs, a laboratory-frame symmetry breaking is observed. To quantify this effect, the experiments aimed at determining the asymmetry parameter, which is formally defined as in equation (61), but with an important difference: at variance with the molecular-frame asymmetries measured in one-photon ionization [71, 86] and discussed in section 4, laboratory-frame asymmetries are determined by considering all electrons ejected in a given direction with respect to the polarization vector, regardless of the proton ejection direction. The dependence of the calculated asymmetry parameter with respect to the proton kinetic energy and the time-delay between the UV and the IR pulse is plotted in figure 14. The time-delay is defined as the elapsed time between the centers of each pulse. Two main mechanisms, labeled mechanism I and II in the figure, were identified as being responsible for the asymmetry. In mechanism I, the asymmetry comes from the interference between autoionization into the lowest molecular ionic state $1s{\sigma }_{g}$ and direct photoionization by the combined action of UV+IR photons. Mechanism I is thus observed for the shortest time delays, when the pulses efficiently overlap in time, favoring the simultaneous absorption of a UV and an IR photon. In mechanism II, there is a laser-induced transition from the $1s{\sigma }_{g}$ state of ${{\rm{H}}}_{2}^{+}$ (resulting from ionization of ${{\rm{H}}}_{2}$ by the UV pulse) to the $2p{\sigma }_{u}$ state of ${{\rm{H}}}_{2}^{+}$. This is a sequential process in which the NWP created by the UV pulse in the ground state of ${{\rm{H}}}_{2}^{+}$ moves until the IR pulse couples it to the first excited ionic state through absorption of a few IR photons. This mechanism can thus be observed even when the pump and probe pulses do not overlap.

Zoom In Zoom Out Reset image size

It is worth stressing that the molecular-frame asymmetry discussed in 4 is strictly due to one-photon absorption processes. It results from the interference between two dissociative paths that involve ionic states of different inversion symmetry (gerade and ungerade) and electrons with different angular momenta, i.e., continuum states of the type {${{\rm{H}}}_{2}^{+}(1s{\sigma }_{g})$^*$\varepsilon {{\ell }}_{\mathrm{odd}}$} and {${{\rm{H}}}_{2}^{+}(2p{\sigma }_{u})$^*$\varepsilon {{\ell }}_{\mathrm{even}}$}. The laboratory-frame asymmetry described here also involves ionic states of difference inversion symmetry, but at variance with molecular-frame asymmetries, they involve electrons with the same angular momentum. This possibility arises from continuum–continuum ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$↔${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ transitions induced by the IR field, which open ionization channels of the form {${{\rm{H}}}_{2}^{+}(1s{\sigma }_{g})$^*$\varepsilon {{\ell }}_{\mathrm{odd}}/{{\ell }}_{\mathrm{even}}$} and {${{\rm{H}}}_{2}^{+}(2p{\sigma }_{u})$^*$\varepsilon {{\ell }}_{\mathrm{odd}}/{{\ell }}_{\mathrm{even}}$} (see [144, 145]).

The proposed mechanisms are confirmed by the NKE distribution as defined in equation (43), i.e., integrated over all the electron ejection angles and energies, and plotted as a function of time delay. Their role has been made evident in more recent calculations using other laser parameters. For instance, the result of an equivalent UV-pump/IR-probe calculation is plotted in figure 15. Specifically, we now use a single attosecond pulse (SAP) (with 400 as of FWHM), centered at 28 eV and intensity $I$ = ${10}^{9}$ W cm⁻² and a 750 nm IR pulse (∼1.65 eV) with duration 9 fs and $I=3\cdot {10}^{12}$ W cm⁻². The calculated total NKE distributions for the dissociative ionization channel are shown in figure 15(a). Panels (b) and (c) show the contributions from each final symmetry: (i) ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ (optically allowed with the absorption of an odd number of photons from the ground state X${}^{1}{{\rm{\Sigma }}}_{g}^{+}$) and (ii) ${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ (with the absorption of an even number of photons from the ground state or an odd number of photons from the previously populated ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ continuum). In panels (d)–(f) and (g)–(i) the contribution of the 1s${\sigma }_{g}$ and 2p${\sigma }_{u}$ ionization thresholds is shown separately. For negative time-delays (the IR reaches the target before the SAP), the distribution is identical to that obtained when only the UV pulse is present, because the number of IR photons required to reach the lowest ionization limit H${}_{2}^{+}(1s{\sigma }_{g})$ from the ground state in the Franck–Condon region is around ten, which is not possible at the IR intensities used in these calculations. For time delays at which the UV and the IR pulses effectively overlap, between approximately −4 to 4 fs, the NKE distribution signal oscillates with half the period of the IR. This is a rather expected behavior, as shown by a number of recent UV+IR experiments, and is the signature of the so-called laser-assisted photoionization [76, 134, 146, 147]. However, in our problem the modulation with the IR is not strictly symmetric with respect to zero time delay because autoionization breaks this symmetry [76]. This is evident in panels (g)–(i), where the NWP moving in the DES favors a sequential process that leads to fragments into the 2p${\sigma }_{u}$ ionic state at NKE in the interval 10–15 eV. These same panels also show the fingerprint of mechanism II discussed above (see figure 14(b)). At time-delays larger than 4 fs, a strong ionization signal appears in the second ionization threshold for low NKE. This reflects the fact that the NWP, initially created by the UV pulse in the $1{\rm{s}}{\sigma }_{g}$ state of ${{\rm{H}}}_{2}^{+}$ and subsequently evolved in time, is strongly coupled to the excited $2{\rm{p}}{\sigma }_{u}$ state of ${{\rm{H}}}_{2}^{+}$ by the IR field when the NWP reaches an internuclear distance of about 3-4 a.u., which is the farthest it can go in 4 fs. This effect was described many years ago and is known as bond-softening [148–150].

Zoom In Zoom Out Reset image size

A different, widely studied problem, in which overlapping attosecond UV and moderately intense fs IR pulses are used, is electron streaking [151–153]. In streaking, the momentum or kinetic energy of the emitted electron is recorded as a function of the time-delay between the UV and the IR pulses. This technique, mostly applied to atomic systems, is commonly used to elucidate the optical characteristics of a SAP from a structureless photoelectron spectra. However, one can use it the other way around to elucidate the target structure from the streaking patterns induced by a known laser pulse. Figure 16 shows the calculated EKE distributions for ${{\rm{H}}}_{2}$ as a function of time delay. For analysis purposes, we superimpose the values of minus the IR vector potential, $-{A}_{{IR}}$(t), in both panels.

Zoom In Zoom Out Reset image size

In the streaking technique, as described for atomic targets, the UV pulse is much shorter than the IR pulse and one can thus consider the photoelectron being emitted at a well-defined time delay. The photoelectron momentum recorded in time follows, within a simple classical approximation, the vector potential of the IR field, ${{\bf{A}}}_{{IR}}$(t) such that, ${\bf{p}}(t)={{\bf{p}}}_{0}-{{\bf{A}}}_{{IR}}(t)$, with ${{\bf{p}}}_{0}$ being the electron momentum at the instant of ejection. Note that in figure 16 we have only included photoelectron emission for an angle 0° with respect to the polarization vector (aligned with the molecular axis). That is to say, we have plotted the result of equation (44) before integrating over a solid angle. Electron emission at angles different from 0° (or 180°) would require to consider geometrical effects that are ignored in a classical picture [154]. The interaction with the remaining ion is also fully neglected in the classical picture, although other extended models can include it to a given extent [155]. Nevertheless, as can be seen in the left panel of figure 16 (labelled 'Total') the probabilities match the expected classical streaking pattern. The left panel shows the total photoionization probability, which in this case includes all dissociative and non-dissociative channels, i.e. it is the probability that results from integrating over all nuclear continuum states and summing over the bound vibrational states associated to each H${}_{2}^{+}$ ionic state. In the right panel of figure 16, we plot instead only the dissociative photoionization contribution, where we observe a rather modified signal due to the larger contribution of autoionization that produces electrons at different times and energies. Again, we observe that probabilities are not symmetric with respect to zero time delay. For negative delays, the probabilities follow the streaking patterns that correspond to electrons moving in a structureless continuum, which is compatible with the fact that only direct ionization is possible. By contrast, for positive delays the signature of autoionization complicates this pattern. Further investigations are needed to explore the possibility of extracting the ionization times from such streaking signal. In fact, studies on electron streaking in molecules are still very scarce and only very recent works seek to explore the effect of a multi-center target in the expected photoionization time-delay [156–158].

We now explore pump–probe schemes in which only UV radiation with relatively low intensity is used. Due to the large value of the corresponding Keldysh parameter, this UV radiation excites and ionizes the ${{\rm{H}}}_{2}$ molecule without perturbing the Coulomb potentials that the particles feel. In this way, the probe step provides information about the intrinsic molecular dynamics of the system and not about the dynamics induced by the probe itself. This is important because in all the UV+IR probe schemes discussed above, most of the observed dynamics is induced by the intense IR pulse, which somewhat hides the field-free intrinsic dynamics of the system, e.g., autoionization.

The first UV pump–UV probe experiments reported in the literature made use of two identical ultrashort pulses of few tens of femtoseconds generated by the free-electron laser in Hamburg (FLASH), and the chosen targets were ${{\rm{H}}}_{2}$ and ${{\rm{D}}}_{2}$ [74, 75]. The sequence of the two UV pulses leads to double ionization of the molecule: the first pulse ionizes ${{\rm{H}}}_{2}$, inducing NWP dynamics in the ground state of ${{\rm{H}}}_{2}^{+}$ and the second pulse maps the NWP motion into the proton kinetic energy spectra resulting from the Coulomb explosion of the molecular ion following the removal of the second electron. In these pioneering experiments, the time resolution was slightly below 10 fs, which was enough to explore the ratio between direct and sequential double ionization processes.

In general, the available x-ray free-electron laser facilities (such as FLASH in Hamburg, LCLS at Stanford, and others) are already able to produce intense XUV pulses, but the shortest pulse durations that they can provide are of the order of a few fs [159], which is still too long to obtain a clean picture of electron dynamics. At present, the most suitable alternative is to use UV or XUV pulses resulting from HHG processes, as explained in the introduction. In fact, the first experimental attempts using XUV-pump/XUV-probe protocols to retrieve electron dynamics in atoms [160] and in the hydrogen molecule [21] were performed recently using HHG techniques that achieve a time resolution of the order of 1 fs. These works measured the total ionization signal as a function of the time-delay between the two pulses, whose Fourier analysis partly reveals the frequencies involved in the pumped target in good agreement with theoretical simulations. Following these works, full dimensional ab initio theoretical calculations performed on He [161] and ${{\rm{H}}}_{2}$ [22, 77] made use of interferometric schemes based on the use of replicas of attosecond UV pulses that allow one to trace the time evolution of electron wave packets. In particular, in a recent theoretical work on ${{\rm{H}}}_{2}$ [77], some of us showed how the amplitudes and phases of the vibronic (vibrational+electronic) wave packet that is created by the pump pulse can be extracted from the total ionization probabilities. The XUV-pump/XUV-probe scheme used in [77] is shown in figure 17, which we discuss briefly below.

Zoom In Zoom Out Reset image size

In this scheme, we use two identical UV pulses of 12.2 eV (or 14 eV) and 2 fs total duration, and intensities of 10¹² W cm⁻². A plot of the twin pulses in the time and frequency domains is shown in figures 3 and 4. We scan the ionization probabilities for time delays from 0 to 100 fs, thus accounting for the longest oscillatory periods of intermediate vibrational states that can be populated by the pump pulse, namely those associated to the electronic singly excited states B and B' ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ in H₂ [21]. As described in section 4, among the possible ionization events we have both dissociative and non-dissociative direct photoionization as well as autoionization from the lowest series of ${{\mathcal{Q}}}_{1}$ DES. In the present UV–UV problem, however, only the ${{\mathcal{Q}}}_{1}$ DES with symmetry ${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ are populated, because this is the only symmetry that can be reached in a two-photon transition. Note that ionization yields are rather large since absorption of the first photon to the B and B' states is a resonant process (see figure 17).

Figure 17 shows that the first pulse can both excite H₂ through a one-photon transition and ionize it through a two-photon transition. More explicitly, the pump pulse produces a molecular wave packet associated to the B and B' ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ states of H₂ and another wave packet in the ${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ electronic continuum. These two wave packets evolve in time until, after a given time delay, the probe comes. The probe pulse can then create a new wave packet in the ${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ continuum following one-photon absorption from the B and B' ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ states. This wave packet will interfere with that generated earlier in the continuum by the pump pulse. In addition, the probe pulse can also create, in the Franck–Condon region, an exact replica of the vibronic wave packets generated by the pump at time delay zero, as shown in the rightmost panel in figure 17. As in the former case, this will lead to interferences between the wave packets created by the pump and probe pulses. The calculated total dissociative photoionization probability as a function of time delay is shown in figure 18 for two different central frequencies of the twin pulses: 12.2 eV (a) and 14 eV (b). As can be seen in figure 17, photons of 12 eV can only excite the B state of H₂, while photons of 14 eV can excite both the B and B'. Thus, it is not surprising that in the second case the interference patterns in the total ionization yields exhibit a more complex structure.

Zoom In Zoom Out Reset image size

For the 12.2 eV pulses in figure 18(a) the ionization yields show a rapidly oscillating pattern superimposed to a slower one. The latter has a maximum at zero time delay (when both pulses overlap in time) and it repeats every 28–30 fs. This time interval is comparable to the averaged vibrational period for the bunch of vibrational states populated by the pump pulse in the B ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ state. In other words, the slow oscillation pattern reveals the nuclear dynamics associated to the molecular wave packet excited by the pump, and it is visible in the ionization yields because of the mentioned interferences introduced by the two-step process. In this sequential process, one photon from the pump pulse populates a given m intermediate vibronic state, and another photon from the probe pulse is absorbed from the m state to the final continuum states. The interferences between different sequential paths, involving different m intermediate states, reveals the $\mathrm{exp}[-i({E}_{m}-{E}_{m^{\prime} })t]$ beatings, which are responsible for the slow oscillations.

For pulses with energy 14 eV in figure 18(b) the wave packet created by the pump pulse is a coherent superposition of vibronic states associated to the B ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ and B' ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ electronic states. Therefore, the calculated spectrum is not that arising from two disentangled NWPs evolving in different electronic states. To extract the different frequency components, we perform the Fourier transform of the dissociative ionization probability functions (for both 12.2 and 14 eV pulses) and plot them in figure 18(c). We observe that the slow beatings described above are associated to phase differences, $\mathrm{exp}[-i({E}_{m}-{E}_{m^{\prime} })t]$, involving vibronic states of the same (or a few neighboring) electronic states, i.e., vibronic states that differ by 0–2 eV. However, the fast beatings seen in figures 18(a) and (b) are due to phase differences that involve both the ground and the B and B' states, $\mathrm{exp}[-i({E}_{m}-{E}_{0})t]$, i.e., energy differences of the order of 11–15 eV. As shown in [77], the complex interference patterns observed in the ionization yield as a function of time delay can be easily interpreted by using the expression:

$$\begin{eqnarray}&&{P}_{f}(\tau )={\left|{a}_{f}\left(1+{e}^{i{\rm{\Delta }}{E}_{f0}\tau }\right)+\displaystyle \sum _{m}{b}_{{fm}}{e}^{i{\rm{\Delta }}{E}_{{fm}}\tau }\right|}^{2},\end{eqnarray} \tag{ 65 }$$

where ${\rm{\Delta }}{E}_{f0}$ = ${E}_{f}-{E}_{0}$ is the energy difference between the ground and final states of energies E₀ and E_f, respectively, ${\rm{\Delta }}{E}_{{fm}}$ = ${E}_{f}-{E}_{m}$ is the energy difference between a given m vibronic (vibrational+electronic) intermediate state of energy E_m and the final continuum state, a_f stands for the two-photon direct ionization amplitude, and b_fm is the amplitude for the sequential process, which includes the excitation from the ground to a given m vibronic state and the time-delayed transition from the m state to the final f state. Equation (65) is exact when the following two premises are fulfilled: (i) the ground state remains with a population close to unity at any time and (ii) the probability of absorbing three or more photons is negligible, so that there are no further continuum–continuum/bound transitions from the final state f.

The square of the second term, which describes the sequential process, gives a $\mathrm{cos}[({E}_{m}-{E}_{m^{\prime} })\tau ]$ term (the slow oscillations shown in both figures 18(a) and (b) corresponding to the low energy peaks in figure 18(c)). For 12.2 eV, the beatings associated to $({E}_{m}-{E}_{{m}^{\prime }})$ only include energy differences between vibrational states associated to the B ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ state, then we obtain signals corresponding to the vibrational periods of this state. For 14 eV, we observe in addition those beatings associated to energy differences between vibrational states within the B' ${}^{1}{{\rm{\Sigma }}}_{u}^{+}$ state, plus the beatings associated to a vibronic state of the B electronic state and a vibronic state of the B' state.

The crossed term between the direct (a_f) and the sequential process (b_fm) resulting from squaring equation (65) leads to a $\mathrm{cos}[({E}_{m}-{E}_{0})\tau ]$ term corresponding to the transition from the ground to the intermediate states. It thus encodes the information about the coupled electron and nuclear dynamics induced by the pump pulse in the latter states (the fast oscillations in the ionization yields). Indeed, the Fourier transform exhibits peaks around 11-15 eV corresponding to each $\mathrm{exp}[-i({E}_{m}-{E}_{0})\tau ]$ electronic phase, with relative intensities that follow the actual values of the amplitudes that describe the transition $0\to m$. This can be seen in figure 18(c) where, along the upper x-axis, we have included the energy positions for all m vibronic states, i.e. each vibrational state associated to the B and B' states.

In conclusion, the Fourier analysis of the total ionization yields measured in a XUV-pump/XUV-probe experiment with twin attosecond pulses allows the full reconstruction of the intrinsic coupled electron–nuclear dynamics induced by the pump pulse. The phases and amplitudes of the pumped molecular wave packet can be easily retrieved from the interferometric effects created by the twin pulses and arising from the different quantum paths that involve two-photon direct and sequential processes. Moreover, as discussed in detail in [77], the different electronic and nuclear timescales involved in the dynamics of the pumped wave packet suggest that, by choosing the appropriate time delay between pump and probe, one could favor the transition from a particular electronic state. This is possible because the nuclear dynamics associated to different electronic states comes with different vibrational periods, so that the resulting coherent molecular wave packet still reconstructs at different times in the inner region, thus leading to a selective localization of the wave packet in a particular electronic state. The main conclusions and mechanisms revealed by the pump–probe schemes discussed here for H₂, as well as most of the analysis performed in the two previous sections, can be easily extrapolated to other diatomic molecules and probably to larger molecular targets, for which full dimensional theoretical tools are not available at present but experiments can be readily performed.