Time delay in XUV/IR photoionization of H$_2$O

TL;DR

This paper presents a computational approach combining coordinate scaling and density functional theory to analyze time delays in XUV/IR photoionization of water molecules, providing insights into ultrafast electron dynamics.

Contribution

It introduces a novel method for solving the time-dependent Schrödinger equation for molecules under combined XUV and IR pulses, specifically applied to water and hydrogen.

Findings

Determined time delays in water and hydrogen photoionization.

Validated the computational approach against experimental data.

Provided detailed temporal electron dynamics in molecular ionization.

Abstract

We solve the time-dependent Schr\"odinger equation describing a water molecule driven by a superposition of the XUV and IR pulses typical for a RABBITT experiment. This solution is obtained by a combination of the time-dependent coordinate scaling and the density functional theory with self-interaction correction. Results of this solution are used to determine the time delay in photoionization of the water and hydrogen molecules.

Tables

Table 1. Table 1: Calculated and experimental ionization potentials of the H 2 O molecule

Orbital	DFT-SIC	Expt.Ning et al. (2008)
$1b_1$	12.3	12.6
$3a_1$	14.9	14.8
$1b_2$	18.3	18.7
$2a_1$	34.3	27.1

Equations

$$i\frac{\partial \Psi (\mathbf{r},t)}{\partial t}=\hat{H}\Psi (\mathbf{r},t)$$

$$\hat{H}=\frac{{\hat{p}}^2}{2}-\mathbf{A}(t)\hat{\mathbf{p}}+U(\mathbf{r}).$$

$$\mathbf{A}(t)=-{\int }_0^{t}q\mathbf{E}(t^{\prime })dt.$$

$$\mathbf{A}(t)={\mathbf{A}}_{\mathrm{XUV}}(t)+{\mathbf{A}}_{\mathrm{IR}}(t-\tau ),$$

$${\mathbf{A}}_{\mathrm{XUV}}(t)=-{\mathbf{n}}_{\mathrm{XUV}}{A}_{\mathrm{XUV}}\exp \left(-2\ln 2\frac{t^2}{{\tau }_{\mathrm{XUV}}^2}\right)\cos {\omega }_{\mathrm{XUV}}t,$$

$${\mathbf{A}}_{\mathrm{IR}}(t)=-{\mathbf{n}}_{\mathrm{IR}}{A}_{\mathrm{IR}}{\cos }^2(\pi t/{\tau }_{\mathrm{IR}})\cos \omega t,|t|<{\tau }_{\mathrm{IR}}/2,$$

$$\Psi (\mathbf{r},t_0)={\phi }_0(\mathbf{r})\exp (-i{E}_0{t}_0),$$

$$\mathbf{r}=a(t)\mathbf{\xi }.$$

$$\Psi (a(t)\mathbf{\xi },t)=\frac{1}{[a(t){]}^{3/2}}\exp \left(\frac{i}{2}a(t)\dot{a}(t){\xi }^2\right)\psi (\mathbf{\xi },t).$$

$$i\frac{\partial \psi (\mathbf{\xi },t)}{\partial t}$$

$$|f(\mathbf{k}){|}^2={\dot{a}}_{\infty }^{-3}\underset{t\to \infty }{\lim }|\psi (\mathbf{k}/{\dot{a}}_{\infty },t){|}^2.$$

$$P^{(3)}\equiv \frac{dP}{d{k}_{x}d{k}_{y}d{k}_z}=|f(\mathbf{k}){|}^2\simeq {\dot{a}}_{\infty }^{-3}|\psi (\mathbf{k}/{\dot{a}}_{\infty },t_f){|}^2.$$

$$\displaystyle a(t)=\left\{\begin{array}[]{ll}1,&t<t_{1};\\ \dot{a}_{\infty}t,&t>t_{1}.\end{array}\right.$$

$$\psi (\mathbf{\xi },t_1+0)=\exp \left(\frac{i}{2}{\dot{a}}_{\infty }{\xi }^2\right)\psi (\mathbf{\xi },t_1-0).$$

$$U_{sa}(\xi ,t)=-i\frac{s(t)}{a(t)}{e}^{-{\xi }^2}$$

$$\displaystyle s(t)=\left\{\begin{array}[]{ll}0,&t<t_{1};\\ (1+\cos[\pi(t-2t_{1})/t_{1}])/2,&t\in[t_{1},2t_{1}];\\ 1,&t>2t_{1};\end{array}\right.$$

$$E_{\text{SIC}}=E_{\text{H}}\{\rho \}+E_{\text{XC}}\{{\rho }_{\uparrow },{\rho }_{\downarrow }\}-[E_{\text{H}}\{{\rho }_i\}+E_{\text{XC}}\{{\rho }_i,0\}],$$

$$\rho (\mathbf{r})=\sum _{i=1}^{N_e}{\rho }_i(\mathbf{r})$$

$${\rho }_i(\mathbf{r})=|{\phi }_i(\mathbf{r}){|}^2.$$

$${\rho }_{\uparrow }(\mathbf{r})={\rho }_{\downarrow }(\mathbf{r})=\rho (\mathbf{r})/2$$

$$E_{\text{H}}\{\rho \}=\frac{1}{2}\int \int \frac{\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })}{|{\mathbf{r}}^{\prime }-\mathbf{r}|}d{\mathbf{r}}^{\prime }d\mathbf{r}$$

$$E_{\text{XC}}\{{\rho }_{\uparrow },{\rho }_{\downarrow }\}=\int [{\rho }_{\uparrow }(\mathbf{r})+{\rho }_{\downarrow }(\mathbf{r})]{\epsilon }_{\text{XC}}[{\rho }_{\uparrow }(\mathbf{r}),{\rho }_{\downarrow }(\mathbf{r})]d\mathbf{r},$$

$$u_i(\mathbf{r})=\frac{\delta E_{\text{SIC}}}{\delta {\rho }_i(\mathbf{r})}$$

$$U(\mathbf{r})=u_{\text{nucl}}(\mathbf{r})+u_i(\mathbf{r})$$

$$-\frac{\partial {\phi }_i(\mathbf{r},t)}{\partial t}=\left[\frac{{\hat{p}}^2}{2}+u_{\text{nucl}}(\mathbf{r})+{\alpha }_k{u}_i(\mathbf{r})\right]{\phi }_i(\mathbf{r},t)$$

$${\phi }_i(\mathbf{r},t+0)={\phi }_i(\mathbf{r},t)-\sum _{j=1}^{i-1}⟨{\phi }_j(\mathbf{r})|{\phi }_i(\mathbf{r},t)⟩{\phi }_j(\mathbf{r})$$

$${\phi }_i(\mathbf{r})=\underset{t\to \infty }{\lim }{\phi }_i(\mathbf{r},t)$$

$${ϵ}_i=-\frac{1}{2}\underset{t\to \infty }{\lim }\left[\frac{1}{⟨{\phi }_i(\mathbf{r},t)|{\phi }_i(\mathbf{r},t)⟩}\frac{d⟨{\phi }_i(\mathbf{r},t)|{\phi }_i(\mathbf{r},t)⟩}{dt}\right]$$

$$S_{2q}(\tau )=A+B\cos [2\omega (\tau -{\tau }_a)],$$

$${\tau }_a={\tau }_{\mathrm{W}}+{\tau }_{\mathrm{CLC}},$$

Full text

I Introduction

Time delay in molecular photoionization can now be measured by using the RABBITT technique (Reconstruction of Attosecond Beating By Interference of Two-photon Transitions). In the pioneering experiment, Huppert et al. (2016) determined the relative delay in photoemission of the outer-most valence shells of the H2O and N2O molecules. The Eisenbud-Wigner-Smith component of the measured time delay (Wigner time delay for brevity), related to the XUV photon absorption, was evaluated from the complex dipole matrix elements provided by molecular quantum scattering theory Baykusheva and Wr̈ner (2017). The Coulomb-laser coupling (CLC) correction, associated with the IR dressing field, was accounted for separately in an atomic-like fashion.

In the present work, we provide an alternative approach where the ionizing XUV field and the dressing IR field drive the molecular time-dependent Schrödinger equation (TDSE). This equation is solved by a combination of the time-dependent coordinate scaling (TDCS) and the density functional theory with self-interaction correction (DFT-SIC). We applied the TDCS technique in our earlier work to describe a RABBITT measurement on the molecular H${}_{2}^{+}$ ion Serov and Kheifets (2016). An advantage of this technique is that the TDSE is solved directly in the XUV and IR fields and the calculated (and measurable) time delay is not artificially split into the Wigner and CLC components. The earlier application to H${}_{2}^{+}$ was based on the explicit form of the one-electron potential. For more complex molecules like H2 and H2O , this potential can be evaluated within the DFT-SIC approach.

II Theoretical methods

II.1 Time-dependent scaling method

We restrict ourselves with a single active electron (SAE) approximation and write the TDSE as

[TABLE]

with the Hamiltonian

[TABLE]

Here $\hat{\mathbf{p}}=-i\nabla$ is the momentum operator, $U(\mathbf{r},t)$ is the electron-nucleus interaction, $\mathbf{A}(t)$ is the vector potential of the electromagnetic field. The latter is defined as 111The atomic units are in use throughout the paper such that $e=m=\hbar=1$. The factor $1/c$ with the speed of light $c\simeq 137$ and the electron charge $q=-1$ are absorbed into the vector potential.

[TABLE]

Here $\mathbf{E}(t)$ is the electric field vector. In a typical attosecond streaking or a RABBITT experiment, the target atom or molecule is exposed to a combination of the two fields:

[TABLE]

where $\tau$ is the relative displacement of the XUV and IR pulses. We model an ultrashort XUV pulse by a Gaussian envelope

[TABLE]

with the FWHM $\tau_{\rm XUV}$. The IR pulse is described by the $\cos^{2}$ envelope

[TABLE]

where $\tau_{\rm IR}$ is the IR pulse duration. The time evolution of the target under consideration starts from the initial state

[TABLE]

where $t_{0}=-\tau_{\rm IR}/2+\tau$ and $\varphi_{0}(\mathbf{r})$, $E_{0}$ are the wave function and the energy of the initial state.

After the end of the XUV pulse, the ionized electron is exposed to a slow varying IR field and the long range Coulomb field of the residual ion. During the propagation in the IR field, the photoelectron gains a considerable speed and travels large distances from the parent ion. To describe this process, solution of the TDSE should be sought in a very large coordinate box for a very long propagation time which places a significant strain on computational resources. To bypass this problem, we employ an expanding coordinate system Serov et al. (2007). In this method, which we term the time-dependent scaling (TDS), the following variable transformation is made:

[TABLE]

Here $a(t)$ is a scaling factor with an asymptotically linear time dependence $a(t\to\infty)=\dot{a}_{\infty}t$ and $\bm{\xi}$ is a coordinate vector. Such a transformation makes the coordinate frame to expand along with the wave packet. In addition, the following transformation is applied to the wave function

[TABLE]

Such a transformation removes a rapidly oscillating phase factor from the wave function in the asymptotic region Serov et al. (2007). Thus transformed wave function satisfies the equation

[TABLE]

where $\hat{\mathbf{p}}_{\xi}=-i\nabla_{\xi}=-i\left(\frac{\partial}{\partial\xi_{x}},\frac{\partial}{\partial\xi_{y}},\frac{\partial}{\partial\xi_{z}}\right)$. Remarkable property of the expanding coordinate system is that the ionization amplitude $f(\mathbf{k})$ is related with the wave function $\psi(\bm{\xi},t)$ by a simple formula Serov et al. (2007)

[TABLE]

In practice, the evolution is traced for a very large time $t_{f}\gg\tau_{\rm IR}$ and then the ionization probability density is obtained from the expression

[TABLE]

The coordinate frame (8) is well suited for approximating an expanding wave packet. However, its drawback is that the bound states are described progressively less accurately as the coordinate frame and its numerical grid expands. Therefore, during the XUV pulse, when an accurate approximation of the bound states is required, we use a stationary coordinate frame. The expansion of the frame starts at the moment $t_{1}\gg\tau_{\rm XUV}$. We use the piecewise linear scaling

[TABLE]

At $t<t_{1}$ the wave function $\psi(\bm{\xi},t)=\Psi(\mathbf{r},t)$. Since the time derivative of $a(t)$ defined by Eq. (15) have discontinuity at the start of the expansion, the wave function at $t_{1}$ should be multiplied by the phase factor

[TABLE]

Here we choose $\dot{a}_{\infty}=1/t_{1}$. Such a choice ensures that the wave packet remains stationary in the expanding frame at $t>t_{1}$.

Bound states are suppressed by introducing an imaginary absorbing potential near the origin:

[TABLE]

Unlike in the previous treatment Serov and Kheifets (2016), here we use a smooth switching of the imaginary absorbing potential, by setting

[TABLE]

This reduces spurious transitions from the bound states to continuum.

After having obtained the ionization amplitude for each individual XUV pulse in the APT, we summed the results. This procedure is correct, because the intensity of the XUV field in RABBITT experiments is usually small, and the processes of sequential absorption of several XUV photons can be neglected.

II.2 Density functional theory with self-interaction correction

We employ the density functional theory (DFT) with the self-interaction correction (SIC) Perdew and Zunger (1981). This correction is necessary to restore the Coulomb asymptotics of the photoelectron interaction with the residual ion which is essential for time delay calculations.

The density functional with the SIC Perdew and Zunger (1981) contains the Hartree $E_{\text{H}}\{\rho\}$ and the exchange-correlation $E_{\text{XC}}\{\rho_{\uparrow},\rho_{\downarrow}\}$ components:

[TABLE]

where the electron density

[TABLE]

is the sum of the particle densities of the $i$-th electron orbital

[TABLE]

Here we consider only molecules with fully coupled electrons and hence the density of the spin-up and spin-down electrons are equal:

[TABLE]

The Hartree energy

[TABLE]

The exchange-correlation functional is expressed in the local density approximation (LDA)

[TABLE]

where $\varepsilon_{\text{XC}}[\rho_{\uparrow}(\mathbf{r}),\rho_{\downarrow}(\mathbf{r})]$ is the exchange-correlation energy per an electron. The effective potential acting upon an $i$-th electron by the rest of the many-electron ensemble is expressed as a functional derivative

[TABLE]

The Kohn-Sham effective potential is the sum of the nuclear and electron components:

[TABLE]

The SIC can be applied to any density functional. In the present application, we used LDA exchange-correlation functional proposed in Gunnarsson and Lundqvist (1976).

Calculation of the one-electron orbitals $\varphi_{i}(\mathbf{r})$ and corresponding one-electro potentials $u_{i}(\mathbf{r})$ is carried out by the imaginary time evolution based on the solution of the equation

[TABLE]

The orthogonality to the occupied states is enforced on each step of the time evolution:

[TABLE]

The stationary orbital is evaluated as the limit

[TABLE]

After finding $\varphi_{i}(\mathbf{r})$, a new set of potentials $u_{i}(\mathbf{r})$ is determined and fed into Eq. (30) to start the next iteration. On the first iteration, the parameter $\alpha_{1}=0$, i.e. only the nuclear term is taken into account in Eq. (30). On the next iterations, $\alpha_{k>1}$ grows linearly with $k$ reaching the value of $\alpha_{k\gg 1}=1$ . This way the inter-electron interaction is switched on gradually thus ensuring a smooth convergence of the solution.

The orbital energies are calculated as

[TABLE]

Thus determined the ionization potential of the H2 molecule is equal to $|\epsilon_{1}|$=16.7 eV whereas the experimental value is 15.6 eV. The calculated ionization potentials $|\epsilon_{i}|$ of the H2O molecule are shown in Table 1 in comparison with the experimental values from Ning et al. (2008).

To solve the TDSE (10) we employed the one-electron effective potential $u_{i}(\mathbf{r})$ corresponding to the ground stationary state and this potential was frozen during the whole ionization process (the frozen-core approximation).

II.3 Numerical implementation

We solve the TDSE. (10) using the orthogonal fast spherical Bessel transform as described in Serov (2017). In all the calculations, we set the box size to $\xi_{max}=51.2$ a.u. The radial grid step was set to $\Delta\xi=0.2$ a.u. The angular basis parameters were $N_{\theta}=16$ and $N_{\phi}=8$ for H2, and $N_{\theta}=4$ and $N_{\phi}=6$ for H2O. The APT is modeled by a series of $N_{\rm APT}=11$ Gaussian pulses with the width $\tau_{\rm XUV}=1$ a.u. (24 as). Such a short pulse duration leads to a large spectral width and allows to obtain the time delay results in a wide range of photoelectron energies.

The APT width $\tau_{\rm APT}=2T_{\rm IR}$ (5.2 fs), whereas a long IR pulse is modeled by a continuous wave with the frequency $\omega=0.05841$ a.u. (photon energy 1.59 eV, $\lambda=780$ nm) and the vector potential amplitude $A_{\rm IR}=0.025$. The amplitude of the XUV pulse was $A_{\rm XUV}=0.2$ a.u. The relative APT/IR time delay $\tau$ was varied from 0 to 0.5$T_{\rm IR}$ with a step 0.03125$T_{\rm IR}$. The central frequency was $\omega_{\rm XUV}=29\omega$ for H2 and $\omega_{\rm XUV}=39\omega$ for water.

By exposing an atom or a molecule to the APT with the central frequency $\omega_{\rm XUV}=(2q_{0}+1)\omega$, the photoelectrons will be emitted with the energies $E_{2q+1}=(2q+1)\omega-E_{0}$ corresponding to the odd harmonics of the IR frequency $\omega$. Superimposing a dressing IR field will add additional peaks in the photoelectron spectrum at $E_{2q}=2q\omega-E_{0}$ These additional peaks, known as the sidebands (SB), correspond to the even harmonics. The sideband amplitudes will vary with the relative time delay $\tau$ of the APT and the IR pulses as Paul et al. (2001)

[TABLE]

where $\tau_{a}$ is the atomic time delay. The atomic time delay can be written in a form

[TABLE]

where $\tau_{\rm W}$ is the Wigner time delay Wigner (1955) and $\tau_{\rm CLC}$ is the Coulomb-laser coupling (CLC) correction Nagele et al. (2011). Here we assume that there is no group delay (chirp) in the APT spectrum and all the harmonics have the same phase.

III Results

In Fig. 1 we display the atomic time delay calculated by the TDCS method. Every energy point on the graph corresponds to a given SB. The error bars indicate the accuracy of the cosine fit to Eq. (34). For comparison, another calculation is shown in which the Wigner time delay is calculated by the prolate spheroidal exterior complex scaling (PSECS) Serov et al. (2013) and the CLC correction is introduced analytically Serov et al. (2015). The PSECS is an ab initio technique and it returns the exact Wigner time delay for diatomic molecules. However, in the photoelectron energy range between 10 and 36 eV, the PSECS calculation on H2 did not converge. Most likely, this is because of a large number of quasi-stationary states in this spectral range.

The TDCS and PSECS+CCLC results agree very well close to the threshold and are qualitatively similar at large excess energies. In the parallel molecular orientation, both set of calculations display a peak in the atomic time delay. However, in the TDCS calculation this peak is shifted by 7 eV towards lower photoelectron energies ($E_{e}=35$ eV in TDCS versus 42 eV in PSECS+CCLC). Such a large difference can be explained by poor performance of the DFT for such a few-electron systems like H2 . We note that the peak displacement by 7 eV far exceeds an error of 1 eV in the ionization potential. This indicates that such a dynamic quantity as the atomic time delay is much more sensitive to inter-electron correlation than the static ionization potential.

Fig. 2 shows the atomic time delay in H2O molecule corresponding to ionization of the HOMO $1b_{1}$ and the HOMO-1 $3a_{1}$. Various plotting symbols display the time delay for randomly oriented and aligned molecules. The randomly oriented results are averaged over all the molecular orientations. The aligned results correspond to photoelectron emission in the direction perpendicular to the nodal plane of the given orbital. In the case of $1b_{1}$ this plane contains the H atoms and the molecular dipole momentum vector. Therefore both perpendicular orientations produce identical results. The nodal plane of $3a_{1}$ is perpendicular to the dipole momentum and hence there are two distinctive perpendicular orientations: in the direction of the H atoms ($\theta=0^{\circ}$) and the reverse direction ($\theta=180^{\circ}$).

In the case of the $1b_{1}$ orbital, the randomly oriented time delay is nearly identical to the aligned results as the perpendicular emission is by far dominant. In the case of the $3a_{1}$ orbital, results are more complicated. At large photoelectron energies, the atomic time delay for $\theta=180^{\circ}$ is large than that for $\theta=0^{\circ}$ which means the escape is faster in the O direction than in the H direction. However, this trend is reversed near threshold. Away from the threshold, the non-aligned time delay is in between these two cases whereas it is above both the aligned curves close to the threshold. This means that some other orientations contribute strongly.

In Fig. 3 we display the time delay difference between the $3a_{1}$ and $1b_{1}$ orbitals at the same photon energy. It may seem surprising that the time delay difference is nearly vanishing in the whole studied energy range. This, however, may be explained by the fact that the randomly oriented water molecule may look like the neon atom and these two states in Ne differ only by the nodal plane orientation. In the same figure we plot the experimental results by Huppert et al. (2016) which show a sign variation of the time delay difference. The frozen-core Hartree–Fock (FCHF) calculation by the same authors Huppert et al. (2016); Baykusheva and Wr̈ner (2017) is also overplotted. Neither calculation reproduce the experimental results within the stated error bars.

IV Conclusions

We combined the time-dependent coordinate scaling (TDCS) method, which we developed earlier for modeling of RABBITT experiments Serov and Kheifets (2016), with the density functional theory with self-interaction correction (DFT-SIC). An advantage of the TDCS method is the direct solution of the time-dependent Schrödinger equation driven by the superposition of the XUV and IR pulses without the splitting of the atomic time delay into the Wigner and CLC components. By using this technique, we calculated the photoionization time delay of the H2 and H2O molecules. In the case of H2 , we made a comparison with the ab initio PSECS calculations and found a good agreement. In the case of H2O , a comparison was made with experiment and agreement was found poor, but other theoretical methods like FCHF did not perform any better.

The theory consistently point to nearly identical time delay from the HOMO and HOMO-1 orbitals of the randomly oriented water molecule. However, this time delay difference becomes noticeable on the oriented H2O molecule. Hence, the experiments will be highly desirable in this oriented configuration which may reveal a reach and anisotropic ultrafast photoelectron dynamics.