The impact of electron–electron correlation in ultrafast attosecond single ionization dynamics

The model we apply in this work represents an extension to the one originally suggested by Shin and Metiu [31, 32]. In their work, a linear, one-dimensional charge-transfer model system was employed consisting of two fixed nuclei with charges Z₁ and Z₂ at the positions ±L/2, one moving nucleus (Z) in between with coordinate R, and one electron, here with coordinate y, giving rise to the potential:

$$\begin{align}\hfill {V}^{1\text{e}}\left(y,R\right)& =\frac{{\mathrm{e}}^{2}}{4\pi {{\epsilon}}_{0}}\left[\frac{{Z}_{1}Z}{\vert L/2-R\vert }+\frac{{Z}_{2}Z}{\vert L/2+R\vert }-\frac{Z\text{erf}\left(\vert R-y\vert /{R}_{\mathrm{c}}\right)}{\vert R-y\vert }\right.\hfill \\ \hfill & \quad -\left.\frac{{Z}_{1}\text{erf}\left(\vert L/2-y\vert /{R}_{\mathrm{f}}\right)}{\vert L/2-y\vert }-\frac{{Z}_{2}\text{erf}\left(\vert L/2+y\vert /{R}_{\mathrm{f}}\right)}{\vert L/2+y\vert }\right],\hfill \end{align} \tag{ 1 }$$

where the error functions (erf) describe a truncated Coulomb interaction between individual particles. The truncation parameters R_f and R_c specify the interaction strength between the electron and the fixed nuclei and the mobile nucleus, respectively [31, 32].

Here, we use an extension to this model introduced by Engel and coworkers, where a second electron, x, is added to the system [33, 34]. The whole particle configuration is shown in figure 1. The system's potential takes on the form

$$\begin{align}\hfill {V}^{2\text{e}}\left(x,y,R\right)& ={V}^{1\text{e}}\left(y,R\right)+\frac{{e}^{2}}{4\pi {{\epsilon}}_{0}}\left[-\frac{Z\text{erf}\left(\vert R-x\vert /{R}_{\mathrm{c}}\right)}{\vert R-x\vert }\right.\hfill \\ \hfill & \quad -\frac{{Z}_{1}\text{erf}\left(\vert L/2-x\vert /{R}_{\mathrm{f}}\right)}{\vert L/2-x\vert }-\frac{{Z}_{2}\text{erf}\left(\vert L/2+x\vert /{R}_{\mathrm{f}}\right)}{\vert L/2+x\vert }\hfill \\ \hfill & \left.\quad +\frac{\text{erf}\left(\vert x-y\vert /{R}_{\mathrm{e}}\right)}{\vert x-y\vert }\right],\hfill \end{align} \tag{ 2 }$$

where R_e scales the electron–electron interaction [34–40]. The fixed nuclei have a distance of L = 10 Å. Nuclear charges are Z = Z₁ = Z₂ = 1. All truncation parameters have been set to R_f = R_c = R_e = 1.5 Å, corresponding to the weak non-adiabatic coupling regime [38], where the electrons adiabatically follow the nuclear dynamics [36]. We want to emphasize that, thus, the model system qualitatively represents a generic molecular environment for the electrons rather than a specific class of molecules, such as linear molecules. Also note that the unscreened interaction of the fixed nuclei with the central nucleus prevents the system from dissociation. Such effects are crucial in the interaction of strong fields with seemingly similar small molecules, such as ${\mathrm{H}}_{3}^{+}$, but are not in the focus of our investigation of excitation and ionization dynamics of molecular systems.

Zoom In Zoom Out Reset image size

The three particle configuration of the model and its dynamics can be fully solved numerically. However, for interpretation, we calculate for the one-electron (1e) and the two-electron (2e) systems the (adiabatic) electronic eigenfunctions, ${\varphi }_{n}^{1\text{e}}\left(y;R\right)$ and ${\varphi }_{n}^{2\text{e}}\left(x,y;R\right)$, with eigenvalues ${V}_{n}^{1\text{e}}\left(R\right)$ and ${V}_{n}^{2\text{e}}\left(R\right)$, respectively, by solving the following eigenvalue equations:

$$\begin{equation}\left[\frac{{p}_{y}^{2}}{2{m}_{\mathrm{e}}}+{V}^{1\text{e}}\left(y,R\right)\right]{\varphi }_{n}^{1\text{e}}\left(y;R\right)={V}_{n}^{1\text{e}}\left(R\right)\;{\varphi }_{n}^{1\text{e}}\left(y;R\right)\end{equation} \tag{ 3 }$$

$$\begin{align}\hfill & \left[\frac{{p}_{x}^{2}}{2{m}_{\mathrm{e}}}+\frac{{p}_{y}^{2}}{2{m}_{\mathrm{e}}}+{V}^{2\text{e}}\left(x,y,R\right)\right]{\varphi }_{n}^{2\text{e}}\left(x,y;R\right)={V}_{n}^{2\text{e}}\left(R\right)\;{\varphi }_{n}^{2\text{e}}\left(x,y;R\right),\hfill \end{align} \tag{ 4 }$$

where m_e is the electron mass and p_x and p_y refer to the electronic momenta. For the two-electron case, the wavefunctions are symmetrized according to Pauli's principle and correspond to an anti-parallel spin configuration. The obtained potential energy curves, ${V}_{n}^{\mathrm{1}\mathrm{e}/\mathrm{2}\mathrm{e}}\left(R\right)$, are shown in figure 2(a) for one (upper panel) and two bound electrons (bottom panel). The one-electron model will be used to analyze the postionization dynamics of one (y) electron remaining in the parent ion after removal of the other (x) electron upon ionization. The first five electronic eigenfunctions of the single electron model, ${\varphi }_{n}^{1\text{e}}\left(y;{R}_{0}\right)$, are shown in figure 2(b) for an asymmetric nuclear configuration, R₀ = −2.05 Å, corresponding to the system initialization (see below). Note, that among these states, for n = 0, 1, and 3, the electron is mostly localized on the left-hand side (y < 0), i.e. at the two close nuclei (in a strongly bound location), whereas for n = 2 and 4 it is predominantly located at the right-hand side (y > 0), i.e. at the single nucleus (weakly bound).

Zoom In Zoom Out Reset image size

The system interacts with an ultrashort attosecond XUV pulse which, using the dipole approximation and velocity gauge, results in the Hamiltonian:

$$\begin{align}\hfill \mathcal{H}\left(t\right)& =\frac{{P}^{2}}{2M}+\frac{{p}_{x}^{2}}{2{m}_{\mathrm{e}}}+\frac{{p}_{y}^{2}}{2{m}_{\mathrm{e}}}+{V}^{2\text{e}}\left(x,y,R\right)\hfill \\ \hfill & \quad +e\;\mathcal{A}\left(t\right)\left(-\frac{Z}{M}P+\frac{{p}_{x}}{{m}_{\mathrm{e}}}+\frac{{p}_{y}}{{m}_{\mathrm{e}}}\right),\hfill \end{align} \tag{ 5 }$$

where M is the proton mass and P the nuclear momentum. The electric field of the ultrashort ionizing XUV pulse is described via its vector potential $\mathcal{A}\left(t\right)$ with a polarization along the molecular axis:

$$\begin{equation}\mathcal{A}\left(t\right)=-\frac{{E}_{0}}{\omega }\;g\left(t-T\right)\;\mathrm{sin}\left(\omega \left(t-T\right)+\phi \right)\end{equation} \tag{ 6 }$$

here, E₀ = −8.7 V Å⁻¹ (or −0.169 au) is the electric field strength (corresponding to an intensity of 1.0 × 10¹⁵ W cm⁻²), ω the field's angular frequency, and g(t − T) a Gaussian pulse envelope function centered around T = 0 fs with a full-width half-maximum (FWHM) of 140 as (2.894 au). The angular frequency corresponds to a wavelength of λ = 80 nm ($\hat{=}\;15.5$ eV or 0.570 au), which is sufficient to singly ionize the molecule through single photon absorption, see figure 2(a). The spectral width of the attosecond pulse intensity is 18.4 eV (0.676 au). The carrier-envelope phase (CEP) ϕ is set to zero, corresponding to a sine-shaped vector potential or an approximately cosine-shaped electric field.

The full system's wave function Ψ(x, y, R, t) is represented on a three-dimensional grid with a range of [−240, 240] Å with 1024 grid points in x- and y-direction, respectively, and of [−4.99, 4.99] Å with 128 points along the R-direction. The time-dependent Schrödinger equation for the Hamiltonian defined in equation (5) is solved numerically with a timestep of 5 as using the split-operator technique [41] and the FFTW 3 library [42] for Fourier transforms. For the details on the numerics please see our previous publications, e.g. reference [38]. The simulation starts at t₀ = −4 fs, well before the XUV pulse enters the system. Reflection at the grid boundaries is suppressed by multiplying Ψ(x, y, R, t) at each timestep with a splitting function [43]

$$\begin{equation}f\left(x,y\right)={\left[1+{\mathrm{e}}^{{\zeta }_{1}\;\left(\vert x\vert -{\zeta }_{2}\right)}\right]}^{-1}\;{\left[1+{\mathrm{e}}^{{\zeta }_{1}\;\left(\vert y\vert -{\zeta }_{2}\right)}\right]}^{-1}\end{equation} \tag{ 7 }$$

with the parameters ζ₁ = 5.67 Å⁻¹ and ζ₂ = 235 Å.

2.2.1. Full fermionic wave function

The initial state is assumed to be a product state of the two-electron adiabatic electronic ground state, ${\varphi }_{0}^{2\text{e}}\left(x,y;R\right)$, and a nuclear wave function

$$\begin{equation}{\Psi}\left(x,y,R,{t}_{0}\right)={\varphi }_{0}^{2\text{e}}\left(x,y;R\right)\;\chi \left(R\right).\end{equation} \tag{ 8 }$$

The adiabatic electronic eigenstates are obtained by solving the field-free electronic Schrödinger equations, equations (3) and (4), via the relaxation method [44].

The nuclear part of the initial wave function, χ(R), is assumed to be a Gaussian-shaped vibrational wave packet centered around the left local minimum of the double-well potential at R₀ = −2.05 Å (see shaded area in figure 2(a)):

$$\begin{equation}\chi \left(R\right)={N}_{0}\;{\mathrm{e}}^{-{\beta }_{R}{\left(R-{R}_{0}\right)}^{2}}\end{equation} \tag{ 9 }$$

here, N₀ serves as normalization constant and the width β_R = 7.14 Å⁻². This way, the Gaussian closely resembles the left-hand side of the vibrational ground state eigenfunction, which is symmetric around R = 0. We note that the factorized, asymmetric initial state does not correspond to the total ground state of the system. It rather corresponds to one of two energetically equal realizations. This is a common situation, where the system resides in one potential well as encountered, for example, in NH₃ inversion or isomerization processes. Restricting the nuclear density asymmetrically to a single potential well gives rise to a specific molecular structure (e.g. one isomer), which allows to distinguish two different binding sites for the electrons, i.e. a strongly bound site (negative electron coordinate) for electrons in the vicinity of the two nuclei at the left-hand side and a weakly bound site (positive electron coordinates) at the single nucleus on the right.

The initial two-electron density, i.e. ∫|Ψ(x, y, R, t₀)|² dR, is displayed in the top layer of figure 3. Note that cuts through the spatial distribution of the electronic part—given by the two-electron groundstate wavefunction ${\varphi }_{0}^{2\text{e}}\left(x,y;R\right)$—approximately correspond to the one-electron functions ${\varphi }_{0}^{\mathrm{1}\mathrm{e}}\left(y;R\right)$ and ${\varphi }_{2}^{\mathrm{1}\mathrm{e}}\left(y;R\right)$. The simulation results obtained for these initial conditions (and their relation to the one-electron electronic eigenfunctions) are discussed in section 3.1. Please note that as we neglect any spin-dependent interaction, spin and spatial coordinates factorize. We therefore only consider the (symmetric) spatial part of the full wave function.

Zoom In Zoom Out Reset image size

2.2.2. Artificial subsystems with distinguishable electrons

For analysis purposes, the full wave function, equation (8), is partitioned into two desymmetrized subsystems with

$$\begin{equation}\left(\mathrm{A}\right)\;\;{\psi }_{\mathrm{A}}\left(x,y,R,{t}_{0}\right)\equiv \sqrt{2}\;{\Psi}\left(x,y,R,{t}_{0}\right)\;{\Theta}\left(y-x\right),\end{equation} \tag{ 10a }$$

$$\begin{equation}\left(\mathrm{B}\right)\;\;{\psi }_{\mathrm{B}}\left(x,y,R,{t}_{0}\right)\equiv \sqrt{2}\;{\Psi}\left(x,y,R,{t}_{0}\right)\;{\Theta}\left(x-y\right),\end{equation} \tag{ 10b }$$

respectively (figure 3, lower panels). In above equation, Θ(x) is the Heaviside step function. These partial wave functions each describing one half of the full system, split along the x = y-diagonal. By doing so, the wave functions ψ_A(x, y, R, t₀) and ψ_B(x, y, R, t₀) vaguely resemble a wave function in Hartree-product form, because now x and y effectively describe identical, yet distinguishable electrons. In the initial configuration of subsystem A (B), the x (y) electron is strongly bound (with an approximate binding energy of ${E}_{B}^{\text{strong}}={V}_{2}^{1\text{e}}\left({R}_{0}\right)-{V}_{0}^{2\text{e}}\left({R}_{0}\right)=15.1$ eV), whereas the y (x) electron (${E}_{B}^{\text{weak}}={V}_{0}^{1\text{e}}\left({R}_{0}\right)-{V}_{0}^{2\text{e}}\left({R}_{0}\right)=10.5$ eV) is weakly bound. Technically, the abrupt cut-off of the wavefunction leads to a weak field-free ionization signal. This background signal is removed from the propagated wave function until t = −250 as, i.e. before the ionizing XUV pulse interacts with the system, by truncating Ψ(x, y, R, t ⩽ −250 as) through multiplication with Θ(25Å − |x|) Θ(25Å − |y|). The different ionization pathways revealed by these subsystems' dynamics are discussed in section 3.2.

2.2.3. Restricted field interaction

For analysis, further disentanglement of individual ionization pathways and their underlying intramolecular dynamics is achieved by artificially restricting the field interaction. To this end, simulations on subsystems A and B are performed with the modified Hamiltonian

$$\begin{align}\hfill {\mathcal{H}}^{\prime }\left(t\right)& =\frac{{P}^{2}}{2M}+\frac{{p}_{x}^{2}}{2{m}_{\mathrm{e}}}+\frac{{p}_{y}^{2}}{2{m}_{\mathrm{e}}}+{V}^{2\text{e}}\left(x,y,R\right)\hfill \\ \hfill & \quad +e\;\mathcal{A}\left(t\right)\left(-\frac{Z}{M}P+\frac{{p}_{\xi }}{{m}_{\mathrm{e}}}\right),\hfill \end{align} \tag{ 11 }$$

where ξ ∈ {x, y} refers to either the x or the y electron. Thus, field interaction is limited to a specific electron. This restriction allows us to distinguish direct and electron correlation driven photoionization pathways. Obtained results are presented and discussed in section 3.3.

To isolate fractions of the wave functions that belong to the singly ionized system, the photoelectron dynamics are evaluated via the ionization signal in regions far from the molecule using the mask function

$$\begin{equation}c\left(x,y\right)={c}_{x}\left(x\right)\cdot \left[1-{c}_{y}\left(y\right)\right]\end{equation} \tag{ 12 }$$

with

$$\begin{equation}{c}_{\xi }\left(\xi \right)=\begin{cases}0\quad \hfill & \mathrm{i}\mathrm{f}\;0\;{< }\;\vert \xi \vert {\leqslant}{\xi }_{\mathrm{c}}\hfill \\ {\mathrm{sin}}^{2}\left(\frac{\left\vert \xi \right\vert -{\xi }_{\mathrm{c}}}{{\Delta}\xi }\frac{\pi }{2}\right)\quad \hfill & \mathrm{i}\mathrm{f}\;{\xi }_{\mathrm{c}}{< }\left\vert \xi \right\vert {\leqslant}{\xi }_{\mathrm{c}}\;+\;{\Delta}\xi \hfill \\ 1\quad \hfill & \mathrm{i}\mathrm{f}\;{\xi }_{\mathrm{c}}\;+\;{\Delta}\xi {< }\left\vert \xi \right\vert {\leqslant}{\xi }_{\text{end}},\hfill \end{cases}\end{equation} \tag{ 13 }$$

where ξ ∈ {x, y} with corresponding x_c = y_c = 25 Å, Δx = Δy = 10 Å, and x_end = y_end = 240 Å marking the end points of the simulation grid. As a result, the mask c(x, y) selects parts of the electronic wave functions

$$\begin{equation}{\psi }_{\text{out}}\left(x,y,R,t\right)\equiv c\left(x,y\right)\;\psi \left(x,y,R,t\right)\end{equation} \tag{ 14 }$$

at large x and low y coordinates, corresponding to the emission of the x electron, whereas the other (y) electron remains bound at the parent molecular ion. Due to symmetry in the electronic coordinates, it is sufficient to only evaluate signals along the x direction. Note, that this mask is independent of R, as the nuclear part of the wave function remains well confined between the two outer nuclei. Further segmentation of ψ_out(x, y, R, t) into subregions S will be introduced in section 3.1.

Interaction of the initial state with the ultrashort XUV pulse induces electron dynamics within the (non-ionized) two-electron system and leads to single ionization. While for atomic systems, an ionization signal can be extracted via projection of the wave function onto a set of Coulomb waves [11], such an approach is not feasible for the multi-centered potential of the molecular model employed here. Instead, we remove the two-electron components for the first 35 two-electron states, for which both electrons are bound, from the total wavefunction. Then we project the remainder (containing only a single bound electron) onto the basis spanned by the electronic eigenfunctions $\left\{{\varphi }_{n}^{1\text{e}}\left(y;R\right)\right\}$ of the one-electron system [2]:

$$\begin{equation}{P}_{n}\left(t\right)=\iint {\left\vert \int {\varphi }_{n}^{1\text{e}}\left(y;R\right)\;{{\Psi}}^{1\text{e}}\left(x,y,R,t\right)\;\mathrm{d}y\right\vert }^{2}\;\mathrm{d}x\;\mathrm{d}R,\end{equation} \tag{ 15 }$$

using the following definitions:

$$\begin{equation}{{\Psi}}^{1\text{e}}\left(x,y,R,t\right)\equiv {\Psi}\left(x,y,R,t\right)-\sum _{m=0}^{34}{a}_{m}^{2\text{e}}\left(R,t\right)\;{\varphi }_{m}^{2\text{e}}\left(x,y;R\right),\end{equation} \tag{ 16 }$$

$$\begin{equation}{a}_{m}^{2\text{e}}\left(R,t\right)\equiv \iint {\varphi }_{m}^{2\text{e}}\left(x,y;R\right)\;{\Psi}\left(x,y,R,t\right)\;\mathrm{d}x\;\mathrm{d}y.\end{equation} \tag{ 17 }$$

The populations P_n(t) are shown in figure 4 upper panel. Note, that in the basis considered here, only contributions are obtained, where the y electron remains bound, while the x electron is ejected. Projection onto $\left\{{\varphi }_{n}^{1\text{e}}\left(x;R\right)\right\}$ yields the same results for reversed roles. One finds that the first five one-electron states are considerably populated through the ultrashort XUV pulse (figure 4 lower panel) around t = 0.

Zoom In Zoom Out Reset image size

For comparison, within the single active electron approximation, the sudden removal of one of the electrons would yield time-independent one-electron state occupations obtained by projection of ${\varphi }_{0}^{2\text{e}}\left(x,y;R\right)$ onto $\left\{{\varphi }_{n}^{1\text{e}}\left(y;R\right)\right\}$:

$$\begin{equation}{b}_{n}^{1\text{e}}\equiv \iint \left\vert \right.\int {\varphi }_{n}^{1\text{e}}\left(y;R\right)\;{\varphi }_{0}^{2\text{e}}\left(x,y;R\right)\;\mathrm{d}y{\left\vert \right.}^{2}\;\mathrm{d}x\;\mathrm{d}R.\end{equation} \tag{ 18 }$$

For the initial conditions considered here, only the states n = 0 (ejecting the weakly bound electron) and n = 2 (ejecting the strongly bound electron) would be populated significantly, leaving the remaining electron more tightly bound in the molecular system (shake-down process [9]). The occupation of higher one-electron states, i.e. n = 1, 3, and 4, corresponding to a shake-up process [9], would be approximately two orders of magnitudes lower (${b}_{1}^{1\text{e}}/{b}_{0}^{1\text{e}}=0.013$, ${b}_{3}^{1\text{e}}/{b}_{0}^{1\text{e}}=0.008$, ${b}_{4}^{1\text{e}}/{b}_{2}^{1\text{e}}=0.006$).

In contrast, in our simulation with fully correlated electrons, depicted in figure 4, these three states show significant occupations. It is also noteworthy, that their population continues to rise after the XUV pulse has passed the system, while in particular the population of the one-electron groundstate (n = 0) declines. Since non-adiabatic transitions occur on a much longer timescale (see below), we trace these phenomena back to the continued interaction between the bound (y) and the ejected (x) electron on early timescales, where a transition to higher bound states is associated with a knock-up process and one to lower states corresponds to a knock-down process [2].

In the following, we aim to disentangle and quantify the various correlation-induced processes that occur during different ionization pathways. To this end, we will evaluate the postionization dynamics. Figure 5 depicts a snapshot of the two-electron density, ∫|Ψ(x, y, R, t)|² dR, of the full antisymmetric system at t = 3 fs after XUV-pulse interaction. As can be gathered, the dominant part of the wave function remains around the origin (corresponding to the non-ionized part of the system), with minor parts being delocalized into the x or y direction. These four double-stripe structures, where either x or y coordinate stays localized within ±10 Å, represent different single ionization processes. Electron densities in regions of high values of both coordinates, |x| and |y|, simultaneously, corresponding to double ionization, are approximately four orders of magnitudes lower due to the much larger energy threshold for double ionization. Consequently, such contributions are not visible in figure 5. In the following, we will concentrate on the single ionization dynamics occurring during the interval of 1 to 20 fs after pulse arrival.

Zoom In Zoom Out Reset image size

Naturally, the electron densities are symmetric with respect to the x = y-diagonal. It is therefore sufficient to restrict the analysis to electron densities emitted along one axis. Here we chose the x axis and ascribe the labels forward/backward for positive/negative values in x. An apparent feature of the single ionization channels is the occurrence of two ionization pathways in every direction. This indicates that the remaining electron eventually stays at different potential minima (around the origin), for which we introduce the labels up/down for positive/negative values in y, respectively, as indicated in figure 5. The electron densities in these four distinct ionization channels differ from each other in shape and amplitude. To distinguish the underlying processes, an evaluation region is defined according to equations (12) and (13) and further separated into subregions S according to the ascribed labels, allowing us to collect and analyze the emitted wave function ${\psi }_{\text{out}}^{S}\left(x,y,R,t\right)$ of each channel separately. As the region of ionization is defined for |x| ⩾ 25 Å, no emission signal is detected until t ≈ 1 fs, when the fastest components of the ionized wave function enter the evaluation region. Ionization signal, i.e. ${I}_{S}\left(t\right)=\iiint {\left\vert {\psi }_{\text{out}}^{\text{S}}\left(x,y,R,t\right)\right\vert }^{2}\;\mathrm{d}x\;\mathrm{d}y\;\mathrm{d}R$, builds up within a subregion S mainly over the period of ∼5 fs. The build up is traced back to a kinetic energy distribution whose central 80% lie between 0.3 and 4.5 eV with a maximum at 1.4 eV. After 5 fs the overall probability to find both electrons in each subregion continues to increase only slightly due to parts of the wave packet with lower kinetic energy entering the subregion. The respective probabilities at 5 fs are 3.9 × 10⁻⁶, 1.1 × 10⁻⁴, 1.8 × 10⁻⁵, and 2.8 × 10⁻⁵ for the forward-up, forward-down, backward-up, and backward-down subregion, respectively.

Figure 6 shows the time evolution of the average momentum, ${\langle {p}_{y}\rangle }^{S}\left(t\right)$, of electron y, which remains bound in the molecular ion after electron x has been emitted. The quantity ${\langle {p}_{y}\rangle }^{S}\left(t\right)$ is calculated via the density distribution for the remaining y electron, ρ_S(p_y, t) by integrating over ${\psi }_{\text{out}}^{S}$ in each individual channels S:

$$\begin{equation}{\langle {p}_{y}\rangle }^{S}\left(t\right)=\frac{\int {p}_{y}\;{\rho }_{S}\left({p}_{y},t\right)\;\mathrm{d}{p}_{y}}{\int {\rho }_{S}\left({p}_{y},t\right)\;\mathrm{d}{p}_{y}},\end{equation} \tag{ 19 }$$

where

$$\begin{equation}{\rho }_{S}\left({p}_{y},t\right)=\iint {\left\vert {\tilde {\psi }}_{\text{out}}^{S}\left({p}_{x},{p}_{y},P,t\right)\right\vert }^{2}\;\mathrm{d}{p}_{x}\;\mathrm{d}P,\end{equation} \tag{ 20 }$$

here, ${\tilde {\psi }}_{\text{out}}^{S}$ is the Fourier-transformed wave function ${\psi }_{\text{out}}^{S}$. The average momentum ${\langle {p}_{y}\rangle }^{S}\left(t\right)$ illustrates that all four regions differ in the dynamics induced in the parent ion.

Zoom In Zoom Out Reset image size

In the following, we will interpret the different dynamics seen in figure 6(a)–(d) based on the leading contributions to photoemission into each evaluation region. However, there are further contributions to each channel, which will be isolated and investigated in sections 3.2 and 3.3.

We interpret the observations as follows: in the forward-down channel, figure 6(d), the signal stems primarily from the direct emission of the weakly bound electron on the right-hand side of the molecule towards positive x values, i.e. without passing the parent ion (and in particular the other electron) first. The strongly bound electron therefore remains on the left-hand side and is hardly affected by the ionization dynamics which is reflected in the nearly constant momentum expectation value of the remaining electron. In contrast, the oscillating signal in the forward-up channel, figure 6(b), can be primarily traced back to the strongly bound electron at negative x values being released towards positive x values, such that it first passes the parent ion and the other electron. Its inelastic scattering with the remaining, weakly bound electron induces oscillations of the latter being reflected by the strong time dependence of the remaining electron's average momentum. We note that the temporal behavior significantly changes for times t > 10 fs, if the nuclear configuration is frozen during the simulation (red dotted lines) suppressing non-adiabatic transitions.

A similar situation but with reversed roles can be seen for emission into the backward direction: here, the down channel, figure 6(c), corresponds to emission of the weakly bound electron after passing the strongly bound one. Thereby, inelastic scattering leads to a regular oscillation in the residual electron's average momentum ⟨p_y⟩(t). Fourier analysis of this oscillation occurring between t = 0 and 10 fs indicates that the corresponding energy of 2.89 eV can be assigned to the energy gap between the electronic ground and first excited state of the one-electron system (2.78 eV at R₀), see figure 2(a). Thus, upon ionization, an electronic wave packet in the residual molecular ion is excited oscillating around the left well's minimum.

The backward-up channel, figure 6(a), on the other hand, corresponds to the direct emission of the strongly bound electron without passing the parent ion. A weak response towards a negative average momentum of the remaining electron, can be seen, in particular between 17 and 20 fs, despite the absence of an immediate interaction between the escaping and the remaining electron. This feature is not present if the simulation is performed with a frozen nuclear configuration (dotted red lines). It originates in the induced nuclear dynamics leading to non-adiabatic transitions (intramolecular charge transfer), which will be discussed in more detail in section 3.3.

The previous analysis provided a first intuitive picture of intramolecular scattering effects in the course of ionization. However, as electrons are indistinguishable, the roles of emitted and remaining electrons during the scattering process cannot clearly be identified. To this end, the artificially truncated wave functions ψ_A and ψ_B, see equations (10a) and (10b), are employed as initial conditions, thus rendering the two electrons distinguishable, see section 2.2.2. This way, electron–electron correlation originating from the antisymmetry of the wave function and interference effects between the two distinct initial density distributions (localized near x = ±5 Å & y = ∓5 Å) are neglected. However, a comparison of the probabilities to find the particles in the evaluation regions, I_S(t), between the full system and the sum of the subsystems A and B shows very good agreement, indicating that for the present system these types of correlation effects are of minor importance.

In figure 7 the time-dependent two-electron densities of the full system with two indistinguishable electrons (upper part) and the subsystems, A (middle part) and B (lower part) with distinguishable electrons, are shown for an area corresponding to the emission of electron x in backward (left panels) and forward direction (right panels), while electron y remains bound to the parent ion. Above the two lower panels, a schematic picture indicates different underlying processes (blue/green arrows) in the x, y-configuration space (left) and the one-dimensional coordinate space (right). The thick blue arrows correspond to the four main contributions, i.e. photoemission with and without intramolecular electron–electron scattering, discussed in section 3.1, where either the strongly (A) or the weakly bound electron (B) interacts with the electromagnetic field and is released to either side of the molecule.

Zoom In Zoom Out Reset image size

The first electron wave packet components enter the evaluation region at x = ±25 Å between 1 and 2 fs after the interaction with the ionizing pulse. At this instant, electron density is mostly found in the backward-up channel for subsystem A and in the forward-down channel for subsystem B (blue arrows) corresponding to direct photoemission from the side of the molecule closest to the respective subregion. Additionally, a strong slightly delayed signal can benoted stemming from an emission into the opposite direction (blue arrows, A: forward-up, B: backward-down), corresponding to photoemission channels involving intramolecular electron–electron scattering, i.e. the x electron first passing the y electron before being finally emitted.

The A/B distinction reveals additional signals with a smaller but still significant probability appearing in the down (A) channels and—to an even smaller extent—in the up (B) channels, which are not visible in figure 5 due to the larger amplitude of the more dominant signals (blue arrows). Their appearance reveals a correlated motion between the two electrons (indicated by green arrows), where the remaining y electron is relocated (intramolecular charge-transfer) either prior to the photoelectron emission of the x electron or afterwards. This question is addressed in the following section 3.3, by further dissecting the ionization pathways through restricting the electrons' interaction with the electric field.

To further investigate the intramolecular dynamics during and after the electron emission process, we perform simulations of the subsystems A and B, i.e. using distinguishable electrons, and restrict the interaction of the electromagnetic field to either the x (ejected) or y (remaining) electron by employing the modified Hamiltonian ${\mathcal{H}}^{\prime }\left(t\right)$ defined in equation (11). This way, the absorption process of a system with distinguishable electrons is strictly limited to one specific single electron.

The ionization wave function, ${\psi }_{\text{out}}^{S}\left(x,y,R,t\right)$, i.e. the part of the wave function entering the analysis region defined via the mask function, equation (13), is projected onto the set of adiabatic eigenfunctions $\left\{{\varphi }_{n}^{1\text{e}}\left(y;R\right)\right\},n\in \left\{0,\dots ,\;4\right\}$ of the one-electron system, obtained from equation (3) and shown in figure 2(b) exemplarily for R = R₀. Note that these states differ in the electron's spatial distributions: while for states with the quantum numbers n ∈ {0, 1, 3}, the y electron is located on the (strongly bound) left-hand side, for n ∈ {2, 4} electron distribution is predominantly found on the (weakly bound) right-hand side. Thus, a transition between states from different sides corresponds to an intramolecular charge transfer. The population ${P}_{n}^{{S}^{\prime }}\left(t\right)$ of the nth one-electron state by the y electron is calculated as

$$\begin{equation}{P}_{n}^{{S}^{\prime }}\left(t\right)=\iint {\left\vert \int {\varphi }_{n}^{1\text{e}}\left(y;R\right)\;{\psi }_{\text{out}}^{{S}^{\prime }}\left(x,y,R,t\right)\;\mathrm{d}y\right\vert }^{2}\;\mathrm{d}x\;\mathrm{d}R,\end{equation} \tag{ 21 }$$

where the domain S' limits the x integration to either positive or negative values corresponding to the forward (S' = fwd) or backward (S' = bwd) channel, respectively, but does not distinguish between up- and down-channels anymore. Thus, ${P}_{n}^{{S}^{\prime }}\left(t\right)$ are the populations of the single-electron states in the evaluation region, i.e. in the molecular parent ion after photoionization. Sketches next to the panels illustrate the dominant ionization pathways with labels '1' and '2' indicating the temporal order.

The populations are shown in figure 8 in backward (left panels, x < 0) and forward direction (right panels, x > 0) for the fully antisymmetric wave function (top panels, a and b) and the subsystems A and B as indicated (lower panels, c–h). The rise in the period of 1 to 5 fs in all panels corresponds to the main part of the wave packet entering the evaluation region. In the subsequent time evolution, the overall probability within the evaluation region remains mostly constant leading to distinct plateau regions in ${P}_{n}^{{S}^{\prime }}\left(t\right)$. However, small contributions with low kinetic energy (of the emitted x electron) continue to enter the evaluation region at later times, while components with high kinetic energy are removed at the grid boundaries. Note that a transient population of metastable states, for example through electronic excitation of 2e states or direct excitation of high-lying vibrational states [45], which further decay into the 1e states, would manifest itself by a delayed ionization signal. Such delayed signals are not observed within the investigated time window. We therefore expect the impact of such states to be negligible in the present case.

Zoom In Zoom Out Reset image size

Several processes (indicated by boldface roman numerals) can be identified and separated from each other. They are summarized with their associated ionization channel S in table 1. We first consider the cases, in which the electron (here, the x electron) interacting with the electric field is the one being eventually emitted (figures 8(c)–(f)). The strongest signal (panel f) corresponds to the emission of the weakly bound electron, initially located near x = +5 Å (subsystem B), in forward direction corresponding to a direct photoemission (I) without scattering with the remaining electron. In this case, the strongly bound y electron remains almost unaffected in its position located at the left-hand side (down), corresponding to the electronic ground state of the one-electron system, n = 0 (see figure 2(b)).

Table 1. Various in parts interdependent processes occur during photoionization and contribute to the signal in different evaluation regions for the two subsystems A (with a strongly bound x electron and weakly bound y electron) and B (vice versa). The build up time (BU) describes the time frame of the initial wave packet entering the analyzing region (defined by signal above the noise level) until a stable population is reached. ${P}_{n}^{{S}^{\prime }}\left(t\right)$ in the last column represents the maximum population of this process.

	Channel S	x	y	Process		Figure 8	BU/fs	${P}_{\text{n}}^{{S}^{\prime }}$/10−5
A	backward-up	<−25 Å	>0	III	Direct emission of the strongly bound electron	c, blue	0.7–5.5	3.1
	backward-down	<−25 Å	<0	VI	Indirect emission via elastic collision with charge	g, green	1.2–14.0	2.5
					transfer and subsequent knock up (inelastic scattering)		1.7–8.2	0.4
				IX	Charge transfer via non-adiabatic transition (following III)	c, green	11.2–19.0	1.0
	forward-up	>+25 Å	>0	IV	Scattered emission of the strongly bound electron	d, blue	1.0–5.9	1.1
					and subsequent knock up (inelastic scattering)		1.3–4.5	0.1
	forward-down	>+25 Å	<0	V	Knock down induced charge transfer (following IV)	d, green	6.2–16.2	0.5
				X	Charge transfer via non-adiabatic transition (following IV)	d, green	11.2–19.0	0.5
	backward-up	<−25 Å	>0	—	Process not visible	—	—	0.0
	backward-down	<−25 Å	<0	II	Scattered emission of the weakly bound electron	e, blue	0.8–7.2	6.6
					and subsequent knock up (inelastic scattering)		1.6–7.0	1.3
B	forward-up	>+25 Å	>0	VII	Indirect emission via elastic collision with charge	h, green	1.5–7.1	0.9
					transfer and subsequent knock up (inelastic scattering)		2.5–4.6	0.04
	forward-down	>+25 Å	<0	I	Direct emission of the weakly bound electron	f, blue	0.7–5.3	22.0
				VIII	Knock-down induced charge transfer (following VII)	h, —	4.2–18.5	0.5
				XI	Charge transfer via non-adiabatic transition (following VII)	h, green	9.9–19.3	0.5

However, the removal of the x electron results in an increase of the y electron's binding energy corresponding to a shake-down process. The energy change corresponding to a sudden electron removal can be estimated from the potential energy curves at R₀, ${\Delta}{E}_{B}\approx {V}^{\text{0e}}\left({R}_{0}\right)+{V}_{0}^{2\text{e}}\left({R}_{0}\right)-{V}_{0}^{1\text{e}}\left({R}_{0}\right)-{V}_{2}^{1\text{e}}\left({R}_{0}\right)$ (where V^0e corresponds to the repulsion energy between the movable and the two fixed nuclei only) and accounts for 1.8 eV. Comparing the relative populations at t = 5 fs, only a very weak shake-up into the one-electron states n = 1 and 3 is noticed with populations of ${P}_{1}^{\text{fwd}}/{P}_{0}^{\text{fwd}}=0.0045$ and ${P}_{3}^{\text{fwd}}/{P}_{0}^{\text{fwd}}=0.0004$ relative to the groundstate. This is in line with the almost constant average momentum seen in figure 6(d) for the full (fermionic) system, indicating no coherent dynamics induced.

In contrast, emitting the weaker bound x electron in the opposite (backward) direction, figure 8(e), such that inelastic scattering with the remaining strongly-bound y electron occurs, entails a significant relative population of the first and third excited one-electron states of ${P}_{1}^{\text{bwd}}/{P}_{0}^{\text{bwd}}=0.21$ and ${P}_{3}^{\text{bwd}}/{P}_{0}^{\text{bwd}}=0.02$, respectively. Note, that these excited states are still localized on the left-hand side of the molecule (down channel). Since such an excitation does not occur in the forward direction, it must be a result of dynamical correlation between the accelerated x electron and the 'inactive' y electron. This interaction corresponds therefore to a pure knock-up process (II) [2]. As a consequence, within 2 and 5 fs, a y electron wave packet can be seen, oscillating within the left potential well, which is reflected in the damped oscillation pattern of the average momentum shown in figure 6(c).

A similar situation is encountered, when the stronger bound electron (figures 8(c) and (d)) interacts with the electric field and is ultimately emitted. If the x electron emission occurs in the backward direction, figure 8(c), again, no intramolecular scattering occurs (III), and the weaker bound y electron remains (initially) in its place, corresponding predominantly to the second excited state, n = 2, which is localized on the molecule's right-hand side (up channel). Again, a shake-down stabilization of the binding energy of approximately 1.8 eV is expected and only a very weak shake-up to state n = 4 is noted (${P}_{4}^{\text{bwd}}/{P}_{2}^{\text{bwd}}=0.0048$).

In the forward direction, i.e. with immediate electron–electron interaction (IV), figure 8(d), the second excited state is dominant, too but also shows a considerable knock-up process (${P}_{4}^{\text{fwd}}/{P}_{2}^{\text{fwd}}=0.10$). Additionally, an increase in the population of the one-electron ground state, n = 0, can be noted. Therefore, inelastic intramolecular scattering with the weaker bound y electron must have taken place resulting in a knock-down process (V) of the y electron during the x electron emission from the (stronger bound) lower energy levels. Comparing their respective peaks, the relative population achieved through the knock-down process is ${P}_{0}^{\text{fwd}}/{P}_{2}^{\text{fwd}}=0.45$. Since in the energetically lower n = 0 state, the y electron is located on the left-hand side, this correlation-driven process coincides with an intramolecular charge transfer (green arrows).

Finally, the lowest panels, figures 8(g) and (h), correspond purely to correlation-driven processes, in which the energy provided by the electric field is absorbed by electron y, but results eventually in the emission of electron x. Therefore, a nearly elastic collision between the two electrons must have occurred, in which the absorbing electron (y) transfers most of its acquired kinetic energy to the electron originally unaffected by the field (x). This 'indirect photoemission' process (VI, VII), is very similar to the two-step-one (TS1) process in double ionization, where the electron emitted first pushes another electron out of an atom in a second step after photoabsorption [10]. But in our case, the initially accelerated electron is not released in the end but rather takes the place of the subsequently emitted electron—similar to the elastic collision between billiard balls. Consequently, the resulting populations of the single-electron states of the remaining electron shown in figure 8(g) for subsystem A resemble the one of subsystem B, seen in panel e, where the remaining electron is located at the strongly bound side. The same applies for subsystem B's elastic collision process, figure 8(h), which rather resembles the populations found in A's inelastic scattering, panel d including knock-up (within VII) and knock-down (VIII) features. Again, this is traced back to the accelerated y electron taking x's place prior to the emission of the x electron. Note that the indirect photoemission leads to a charge transfer (through elastic collision) immediately after photoabsorption. We therefore conclude that the significant early signals, seen in figure 7 in subsystem A's backward-down channel (and with a lower amplitude also in subsystem B's forward-up channel) correspond to the elastic collision process preceding the x electron's emission.

In panels c, d, and h of figure 8, where the n = 2 state is predominantly occupied, a decrease occurs in ${P}_{2}^{{S}^{\prime }}\left(t\right)$ after approximately t ≈ 15 fs with a simultaneous increase of the n = 1 state, ${P}_{1}^{{S}^{\prime }}\left(t\right)$. This time-dependent feature can be traced back to the non-adiabatic nuclear reorganization dynamics (IX, X, XI) induced by the ionization process. Note that these transitions do not occur, if the simulation is performed with a frozen nuclear configuration (dotted lines). Figure 9 shows the correlated electron–nuclear dynamics for subsystem A with emission in backward direction (corresponding to figure 8(c)) through the density function $\int {\left\vert {\psi }_{\text{out}}^{\text{bwd}}\left(x,y,R,t\right)\right\vert }^{2}\;\mathrm{d}x$. It can be seen that during the first 14 fs, the shape of the electronic part only marginally changes, while the center of the nuclear distribution moves from R₀ = −2.05 Å towards larger values. This dynamics is induced by the Coulomb attraction between the remaining electron and the mobile nucleus, but is also consistent with the potential energy surface of the second state, ${V}_{2}^{1\text{e}}\left(R\right)$, see figure 2(a) upper panel. The latter one exhibits a large gradient towards the molecular center (R = 0), where a coupling region with the first excited state, ${V}_{1}^{1\text{e}}\left(R\right)$, is found. Indeed, at 17 fs, the center of the nuclear distribution passes the origin and the electronic distribution begins to shift towards the left-hand side (to negative y values), which is reflected in a slightly negative instantaneous average electronic momentum ${\langle {p}_{y}\rangle }^{\text{backward}-\text{up}}\left(t\right)$ in the case of direct photoionization in the backward-up channel, cf figure 6(a), which is not present in the case of frozen nuclei (red dotted line). Therefore, upon ionization, nuclear dynamics is initiated, driving the system via non-adiabatic transition from the second to the first electronically excited state, ${\varphi }_{2}^{1\text{e}}\left(y;R\right)\to {\varphi }_{1}^{1\text{e}}\left(y;R\right)$ corresponding to an intramolecular charge transfer. Note, however, that this process is significantly slower than the charge transfer process driven by electronic correlation.

Zoom In Zoom Out Reset image size

It is noteworthy that, here, despite these non-adiabatic transitions, the nuclear dynamics does not seem to affect the various ionization processes discussed above. We attribute this to the chosen near-equilibrium initial conditions for χ(R). Previous studies have shown, that nuclear dynamics following initial non-equilibrium configurations is reflected in the photoelectron momentum distribution [37, 38]. Investigations of such effects on the correlation-driven knock-up and knock-down processes are currently under way in our workgroup.

To summarize, all observed processes, i.e. direct photoemission with and without inelastic scattering leading to knock-up and knock-down transitions, as well as indirect photoemission through elastic collision, and non-adiabatic transitions, are summarized in table 1 together with their individual amplitudes at their respective maximum and build up times, i.e. the time span the corresponding wave packet requires to achieve approximately stable populations within the observed time window. These times differ for the individual processes mostly due to the different traveling distances (direct emission vs. scattered and indirect emissions) for the ejected electron to reach the evaluation zone and because of differences in their kinetic energies. However, the slightly longer timescales of knock-down processes in particular indicate a more complex electron–electron dynamics within the molecular system prior to the electron release.

Regarding the relative amplitudes, we note, that the correlation-driven pathways through elastic and inelastic scattering appear to be nearly of the same order of magnitude as the direct photoemission. This can be seen in figures 8(a) and (b), where the dynamics of the fully antisymmetric initial wave function without any restrictions on the electric–field interaction (solid lines) is qualitatively reproduced by the artificially restricted subsystems (dashed lines, corresponding to the direct sum of all individual pathways with restricted field interaction). Remaining discrepancies stem most likely from the missing simultaneous interaction of the XUV pulse with both electrons and also from the omission of interference effects between emitted density from the two initial localized electronic density distributions, which we have dropped by regarding electrons as distinguishable particles.

We conclude that the final state of the molecular parent ion after ionization depends strongly on dynamical correlation between electrons through elastic and inelastic intramolecular scattering events preceding the electron emission beyond static shake-up effects. Furthermore, we showed that such effects also significantly contribute to the total ionization probability. In particular, quantifying the effect of correlation-driven knock-up/knock-down processes and the indirect photoemission processes on the same magnitude as direct photoionization processes, underlines the deficiencies of the commonly applied single active electron approximation as well as the sudden approximation—even in the context of single-photon ionization of molecules.