Ionisation of H$_2$O by a strong ultrashort XUV pulse: a model within the single active electron approximation

TL;DR

This paper introduces a new computational method within the single active electron approximation to study ionisation of water molecules by intense ultrashort XUV pulses, providing efficient and comparative analysis of molecular orientation and phase effects.

Contribution

A novel, computationally inexpensive approach for modeling ionisation in complex molecules under ultrashort laser pulses within the single active electron framework.

Findings

Ionisation yield varies with molecular orientation and pulse phase.

The new method produces results comparable to existing approaches.

Application to water molecule demonstrates practical utility.

Abstract

We present and discuss a new computationally inexpensive method to study, within the single active electron approximation, the interaction of a complex system with an intense ultrashort laser pulse. As a first application, we consider the one photon single ionisation of the highest occupied molecular orbital of the water molecule by a laser pulse. The ionisation yield is calculated for different orientations of the molecule with respect to the field polarization axis and for different carrier envelope phases of the pulse, and compared against predictions of another single active electron approach.

Equations

$$\mathbf{A}(t)=A(t)\mathbf{e}=A_0{\sin }^2\left[\frac{\pi t}{T}\right]\sin (\omega t+\varphi )\mathbf{e},0\le t\le T,$$

$$A_0=\frac{1}{\omega }\sqrt{\frac{I}{I_0}},$$

$$b(t)=-\frac{1}{c}{\int }_0^{t}A(\xi )\mathrm{d}\xi ,$$

$$\left[\mathrm{i}\frac{\partial }{\partial t}-\frac{p^2}{2}-\frac{1}{c}A(t)(\mathbf{e}\cdot \mathbf{p})-\frac{1}{2c^2}{A}^2(t)\right]\Psi (\mathbf{p},t)$$

$$-\int \frac{\mathrm{d}\mathbf{u}}{(2\pi {)}^3}V(\mathbf{p}-\mathbf{u})\Psi (\mathbf{u},t)=0.$$

$$V(\mathbf{p}-\mathbf{u})=-\frac{4\pi }{|\mathbf{p}-\mathbf{u}{|}^2}\approx V(\mathbf{p},\mathbf{u})=-\sum _{n=1}^N{v}_n(\mathbf{p})v_n^{\ast }(\mathbf{u}).$$

$$\Psi (\mathbf{p},t)=e^{-i\zeta (t)}\Phi (\mathbf{p},t),$$

$$\zeta (t)=\frac{1}{2c^2}{\int }_0^{t}A(\xi {)}^2\mathrm{d}\xi .$$

$$F_n(t)=\int \frac{\mathrm{d}\mathbf{u}}{(2\pi {)}^3}{v}_n^{\ast }(\mathbf{u})\Phi (\mathbf{u},t).$$

$$\Phi (\mathbf{p},t)=e^{-\mathrm{i}\frac{p^2}{2}t+\mathrm{i}b(t)p_z}[\Phi (\mathbf{p},0)+\mathrm{i}\sum _{n=1}^N{v}_n(\mathbf{p})$$

$$\times {\int }_0^t\mathrm{d}\xi F_n(\xi )e^{\mathrm{i}\frac{p^2}{2}\xi -\mathrm{i}b(\xi )p_z}].$$

$$\mathbf{F}(t)={\mathbf{F}}_0(t)+{\int }_0^t\mathbf{K}(t,\xi )\mathbf{F}(\xi )\mathrm{d}\xi .$$

$${\Phi }_{\alpha }(\mathbf{r})=\sum _{j=1}^3\sum _{{\gamma }_j}{C}_{{\gamma }_j,\alpha }{\mathcal{G}}_{{\gamma }_j}(\mathbf{r}-{\mathbf{R}}_j),$$

$$V(\mathbf{p},\mathbf{u})=-v(\mathbf{p})v^{\ast }(\mathbf{u})$$

$$\left[{\epsilon }_0-\frac{p^2}{2}\right]{\Phi }_0(\mathbf{p})-\int \frac{\mathrm{d}\mathbf{u}}{(2\pi {)}^3}V(\mathbf{p},\mathbf{u}){\Phi }_0(\mathbf{u})=0,$$

$$a=\int \frac{\mathrm{d}\mathbf{u}}{(2\pi {)}^3}{v}^{\ast }(\mathbf{u}){\Phi }_0(\mathbf{u}),$$

$$v(\mathbf{p})=-\frac{1}{a}\left[{\epsilon }_0-\frac{p^2}{2}\right]{\Phi }_0(\mathbf{p}).$$

$$\left[\mathrm{i}\frac{\partial }{\partial t}-\frac{p^2}{2}+\frac{A(t)}{c}(\mathbf{e}\cdot \mathbf{p})\right]\Phi (\mathbf{p},t)-\int \frac{\mathrm{d}\mathbf{u}}{(2\pi {)}^3}V(\mathbf{p},\mathbf{u})\Phi (\mathbf{u},t)=0,$$

$$\Phi (\mathbf{p},0)={\Phi }_0(\mathbf{p}).$$

$$\Phi (\mathbf{p},t)=e^{-\mathrm{i}{p}^2/2t+\mathrm{i}b(t)(\mathbf{e}\cdot \mathbf{p})}\left[{\Phi }_0(\mathbf{p})+\mathrm{i}v(\mathbf{p}){\int }_0^{t}F(\xi )e^{\mathrm{i}{p}^2/2\xi -\mathrm{i}b(\xi )(\mathbf{e}\cdot \mathbf{p})}\mathrm{d}\xi \right],$$

$$F(t)=\int \frac{\mathrm{d}\mathbf{u}}{(2\pi {)}^3}{v}^{\ast }(\mathbf{u})\Phi (\mathbf{u},t).$$

$$F(t)=F_0(t)+{\int }_0^{t}K(t,\xi )F(\xi )\mathrm{d}\xi .$$

$$S(x,y)$$

$$T(x,y)$$

$$Q(x,y)$$

$$F_0(t)=-\frac{1}{a}{\left[{\epsilon }_{0}S(x,y)-\frac{1}{2}T(x,y)\right]|}_{x=t/2,y=b(t)},$$

$$K(t,\xi )=\frac{\mathrm{i}}{a^2}\left[{\epsilon }_0^{2}S(x,y)-{\epsilon }_{0}T(x,y)+\frac{1}{4}Q(x,y)\right]{|}_{x=(t-\xi )/2,y=b(t)-b(\xi )}.$$

Full text

UCL]Institute of Condensed Matter and Nanosciences, Université Catholique de Louvain, 2 chemin du cyclotron, BOX L7.01.07, B1348 Louvain-la-Neuve, Belgium

\alsoaffiliationSkobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University, Moscow, Russia

UCL]Institute of Condensed Matter and Nanosciences, Université Catholique de Louvain, 2 chemin du cyclotron, BOX L7.01.07, B1348 Louvain-la-Neuve, Belgium

RHUL]Department of Mathematics, Royal Holloway, University of London, Egham, Surrey TW20 0EX, United Kingdom

RHUL]Department of Mathematics, Royal Holloway, University of London, Egham, Surrey TW20 0EX, United Kingdom

UTR]Dipartimento di Scienze Chimiche e Farmaceutiche, Universita’ di Trieste, Trieste, Italy

UTR]Dipartimento di Scienze Chimiche e Farmaceutiche, Universita’ di Trieste, Trieste, Italy

UCL]Institute of Condensed Matter and Nanosciences, Université Catholique de Louvain, 2 chemin du cyclotron, BOX L7.01.07, B1348 Louvain-la-Neuve, Belgium

\abbreviationsSAE,TDSE,DFT

1 Introduction

The recent development of coherent sources of light such as attosecond lasers 1, high-order-harmonic generation (HOHG) sources 2, 3 or free-electron lasers 4, 5 has opened the route to the study of the interaction of matter with intense femtosecond and even sub femtosecond radiation pulses in the XUV regime. Such studies allow one to analyze electron dynamics in atoms or electron and nuclei dynamics in molecules with an unprecedented degree of temporal and spatial resolution. Within this context, the interaction of the water molecule with such radiation pulses is of particular relevance in medical physics, in particular, radiotherapy, since water is one of the main components of most living tissues.

In the present contribution, we develop a new computationally inexpensive method to study, within the single active electron (SAE) model and in the non-relativistic regime, the single ionisation of atoms and molecules by an intense femtosecond or sub femtosecond XUV pulse. Specifically we aim to address the total ionisation yield as function of field frequency, intensity and orientation dependence of the polarisation vector with respect to the molecule. We apply this approach to the water molecule in its ground state while paying attention to the validity of the assumptions we make in such a treatment. We assume that the pulse duration is small compared to the characteristic time of the vibrational motion of the molecule and that, during the interaction, the geometry of the molecule is not modified as we employ the fixed nuclei approximation.

Recently, by measuring the ratio of H2O and D2O high harmonic yields, Farrell et al. 6 managed to characterize the nuclear motion in the molecular states of H2O+. They showed that by contrast to the ionisation of the highest occupied molecular orbital (HOMO), the single ionisation of the second least bound orbital HOMO-1 triggers a fast nuclear dynamics of the molecular ion through a strong bending motion of the molecule. As a result, we only consider here, the single ionisation of the orbital HOMO 1b1 which leaves the geometry of the molecule practically unchanged during the interaction with the pulse. In fact, the period of the fastest oscillation in the water molecule, namely the asymmetrical stretching of the OH bonds, is 8.9 fs 7 which is much longer than the pulse durations we consider here. In other words, we can assume that the molecule is frozen during its interaction with the pulse. However, it is important to note that experimentally, it is impossible to know, a priori, from which orbital the electron is ejected. Farrell’s results show that ionisation from HOMO-1 3a1, which sends the molecule into the $\tilde{A}$ ${}^{2}A_{1}$ state of H${}_{2}O^{+}$, strongly excites the bending mode at photon energies around 14.8 eV. This puts some limitation on the photon energy used in the present work and requires us to pay attention to the bandwidth of the pulse.

Our approach is based on a model that was first developed to treat the interaction of atomic hydrogen with an electromagnetic pulse 8. We work in the momentum space and use the velocity gauge. The main idea is to replace the kernel of the Coulomb potential by a sum of $N$ symmetric separable potentials, each of them supporting a bound state of the system. This method, which we denote by SPAM for Separable Potentials for Atoms and Molecules, allows one to reduce the 4-dimensional time-dependent Schrödinger equation (TDSE) to a system of $N$ 1-dimensional Volterra integral equations depending only on time. As a result, the integration over the spatial coordinates which, in some cases, requires prohibitively large grids or bases, is completely avoided. Each separable potential may be calculated from the exact wave function of the atomic state it supports. However, its analytical expression is not always unique. We developed a procedure to calculate these separable potentials. It provides results for the electron energy spectra that compare very well with those obtained by solving the TDSE with the exact Coulomb potential in situations where the number of essential atomic states playing a significant role is low. By moving from the momentum space to the configuration space, it is easy to show that the separable potentials have a finite range. Let us note that once the separable potentials are determined, the continuum states are automatically defined and, being solutions of the same equation as the one satisfied by the exact bound states taken into account in the calculations, they are orthogonal to these bound states. To generalize to more complex systems such as the water molecule, we proceed along the same lines. We first generate the HOMO in terms of gaussian type orbitals by means of the well established quantum chemistry software package GAMESS(US). It is then straightforward to move to the momentum space and to define the corresponding separable potential which is unique in this case. As for atomic hydrogen, the final step involves solving a 1-dimensional Volterra integral equation. For the sake of completeness, it is important to mention that our approach is not gauge invariant as it is the case for most of the approximate treatments. The problem of the gauge and the delicate question of the possible existence of a privileged gauge are discussed in detail in the context of the present model by Galstyan et al. 9.

The problem of the interaction of the water molecule with a femtosecond or sub femtosecond XUV pulse has almost never been treated theoretically up to now. As far as we know, most of the theoretical calculations have been performed at a wave length of 800 nm. However, it is interesting from the methodological point of view to mention three different contributions. In the first one, Borbély et al. 10 study the ionisation of the water molecule by an intense ultrashort half-cycle electric pulse. They performed both quantum mechanical and classical calculations within the single active electron and frozen core approximation. They considered the ejection from the HOMO 1b1. Since this orbital is mainly constructed from the $2p_{z}$ orbital of the oxygen atom when the molecule is lying in x-y plane, they modelled it by an hydrogenic $2p_{z}$ orbital with an effective charge chosen to reproduce the ionization energy of the HOMO. They found good agreement between the classical and quantum mechanical calculations at high field intensity where the over-the-barrier ionisation mechanism is dominant. In the second contribution, Della Picca et al. 11 study the orientation-dependent single ionisation of fixed-in-space H2O by a short laser pulse for two wave lengths: 800 nm and 76 nm well in the XUV regime. In their calculation, which is based on the strong field approximation (SFA)12 , the initial and final states are described by single-Slater determinants of spin-orbitals, the spatial part of it being calculated by means of the same quantum chemistry software package as in our case. They take into account the five occupied molecular orbitals. The SFA is a first order theory in the sense that the ionisation results from the absorption of a "virtual" photon that is supposed to describe tunneling emission. In the high frequency regime, we have shown 12 that the SFA gives qualitatively good results. In their contribution, Della Picca et al. showed that the HOMO-1 dominates the single-electron emission process when the laser is polarized along the symmetry axis of the molecule and that the electron emission is in general favored in the direction along the laser polarization direction. In the third contribution Petretti et al. 13 apply the single active electron approximation time dependent Schrödinger equation (SAE-TDSE) method to the water molecule. They solve a 4-dimensional TDSE within the single-active-electron approximation to treat the orientation-dependent ionisation of H2O in few-cycle 800 nm linear-polarized laser pulses. The molecular orbitals are Kohn-Sham orbitals obtained by using the LB94 exchange-correlation functional. They showed that the HOMO dominates the overall ionsation behaviour except in the nodal plane of this orbital where the dominant contribution comes from the HOMO-1. The role of the carrier envelope phase is also investigated.

Our contribution is organized as follows. After this introduction, we present our method. First, very briefly in the case of atomic hydrogen and then, in more detail, in the case of the water molecule. In the third section, we first present some tests of our method in the case of atomic hydrogen. Before the conclusions and perspectives we present our results for the water molecule and compare them against the predictions of SAE-TDSE method described in detail by Petretti et al13.

Atomic units combined with the Gaussian system for the electromagnetic field are used throughout unless otherwise specified.

2 Theoretical model

In this section, we define the pulse and describe our model. For the sake of clarity, we first consider briefly the case of atomic hydrogen. Details of the calculations are found in Galstyan et al. 9. We then show, in more detail, how it can be generalized to a more complex system such as the water molecule.

2.1 Description of the pulse

We use the dipole approximation and assume that the electric field is linearly polarized along the unit vector $\mathbf{e}$, the direction of which, unless explicitly stated, coincides with our z-axis. The case of an arbitrary polarization does not introduce any additional complication and does not require more computer resources. The electric field is defined in terms of the vector potential $\mathbf{A}(t)$ as follows:

[TABLE]

where $\omega$ is the radiation frequency, $\phi$, the carrier envelope phase and $T$, the full pulse duration. In term of the peak intensity, the amplitude $A_{0}$ is given by:

[TABLE]

where $I_{0}=3.5\times 10^{16}$ W/cm2 is the atomic unit of intensity. Before considering the atomic hydrogen case, it is convenient to define the following quantity:

[TABLE]

where $c$ is the speed of light.

2.2 Atomic hydrogen

We work in the momentum space and use the velocity gauge. Under these conditions, the TDSE that describes the interaction of atomic hydrogen with the radiation pulse is:

[TABLE]

[TABLE]

$\mathbf{p}$ is the canonical momentum and $V(\mathbf{p}-\mathbf{u})$ is the kernel of the Coulomb potential. The main idea of the model is to replace this kernel by a sum of $N$ symmetric separable potentials supporting $N$ bound states of atomic hydrogen:

[TABLE]

In order to solve Eq. (4) with the kernel of the Coulomb potential given by Eq. (5), we first perform a contact transformation of the wave function $\Psi(\mathbf{p},t)$ to eliminate the $A^{2}(t)$ term from Eq. (4),

[TABLE]

where,

[TABLE]

We then define the following function:

[TABLE]

Under these conditions, the solution of Eq. (4) with the kernel of the Coulomb potential replaced by the expression (5) may be formally written as follows:

[TABLE]

[TABLE]

Upon substituting expression (9) of $\Phi(\mathbf{p},t)$ into Eq. (8), we obtain a system of $N$ coupled linear 1-dimensional Volterra integral equations of the second kind which may be written in matrix form as:

[TABLE]

$\mathbf{F}(t)$ is a vector of dimension $N$, the components of which are the $F_{n}(t)$ functions. $\mathbf{F}_{0}(t)$ is a vector of the same dimension which results from the contribution of the initial wave function $\Phi(\mathbf{p},0)$ present in Eq. (9). $\mathbf{K}(t,\xi)$ which is the kernel of the Volterra equation, is a $N\times N$ matrix (see 8 for the details of the calculations). Consequently, the 4-dimensional TDSE is reduced to a system of $N$ coupled 1-dimensional Volterra integral equations depending only on time which is a major advantage when the integration of the TDSE over the spatial coordinates requires a large grid or basis set.

2.3 Water molecule

To generalize the previous model to the case of the water molecule, we first generate the spatial part $\Phi_{0}(\mathbf{r})$ of the HOMO 1b1 in the configuration space. It is obtained by geometry optimization with the GAMESS(US) program 14 in the Hartree-Fock approximation. Instead of Hartree-Fock one could use DFT, it would not change much, as long as one uses experimental orbital energies instead of ones generated by GAMESS(US) as they are strongly influenced by the method while the wave functions are not. The general expression of a molecular orbital $\alpha$ denoted by $\Phi_{\alpha}(\mathbf{r})$ is:

[TABLE]

where index $j$ designates each nucleus in the molecule. The second sum is over the atomic orbitals taken into account around each nucleus and $\mathcal{G}_{\gamma_{j}}$ is a so-called contracted gaussian from a 6-31G basis set 15 in the present case. $\mathbf{r}-\mathbf{R}_{j}$ is the electronic coordinate relative to the nucleus $j$. The coefficients $C_{\gamma_{j},\alpha}$ are the ones generated by the GAMESS(US) program. The HOMO molecular orbital $\Phi_{0}(\mathbf{r})$ is then expressed in momentum space by Fourier transform. The corresponding separable potential

[TABLE]

is calculated by imposing that the molecular orbital $\Phi_{0}(\mathbf{r})$ is a solution of the following stationary Schrödinger equation:

[TABLE]

where $\varepsilon_{0}=-0.463$ a.u. is the experimental value of the energy of the HOMO 1b1 16. If we define the coefficient

[TABLE]

we obtain immediately:

[TABLE]

The coefficient $a$ is easily obtained by multiplying Eq. (13) by $\Phi_{0}^{*}(\mathbf{p})$ and integrating over $\mathbf{p}$. Consequently, the TDSE that describes, within the SAE approximation, the ejection of an electron from the HOMO 1b1 of the water molecule exposed to a laser pulse is:

[TABLE]

with the initial condition:

[TABLE]

As in the case of atomic hydrogen, the solution of this equation is obtained analytically:

[TABLE]

and the function $F(t)$ is defined by:

[TABLE]

This function is the solution of the following 1-dimensional Volterra integral equation:

[TABLE]

As we stressed before, this function depends only on time. All the spatial dependence is "hidden" in the kernel $K(t,\xi)$. In other words, the spatial dependence is treated analytically which allows one to avoid all the problems related to the size of the grid or of the basis used to integrate on the spatial coordinate of the electron. $F_{0}(t)$ and $K(t,\xi)$ can be expressed in terms of the following functions (see the Appendix for more details):

[TABLE]

In the case of the HOMO 1b1, we have:

[TABLE]

and

[TABLE]

The main shortcomings of the present method are the fact that the separable potential, which is unique if only one molecular orbital of a given symmetry is treated, does not necessarily have the correct asymptotic behaviour as well as the gauge invariance problem. The problem of the gauge invariance in the context of nonlocal potentials has been discussed in great detail in various references 17, 18, 9. In particular, it is shown in 9 that it is possible to formulate our separable potential model in various ways that can be grouped into two families such that within a given family, the length and velocity gauge formulations give the same value for the observables. This shows clearly the non-existence of a privileged gauge but does not solve the problem of the "global" gauge invariance since formulations belonging to the first and the second family give different results for the observables. In the case of atomic hydrogen exposed to a relatively high frequency (of the order or higher than the ionisation potential) field, the present velocity gauge formulation of our model with several bound states taken into account gives results that are in good agreement with the TDSE results with the full Coulomb potential.

2.4 Numerical implementation

In SPAM, Eq. (20) is solved using the block-by-block method 19. From the computational point of view, the most expensive step is the calculation of the kernels $K(t,\xi)$, and it strongly depends on the number of gaussians in the basis set. This number depends on the system under consideration and on the basis set. Calculation for atomic hydrogen with a huge basis set can be as expensive as a calculation for a water molecule, but with a small basis set.

However, the algorithm can be effectively parallelised. We use NVIDIA CUDA for General Purpose Graphics Processor Units (GPGPUs) combined with Message Passing Interface (MPI), so we use several nodes with GPUs. The overall time to calculate the ionisation yield at a high frequency for a water molecule using the 6-31G or 6-311++G** basis set 15 is around 30 minutes and 6 hours respectively on one GPU of CUDA compute capability 2.0. We use two GPUs with MPI to control the error during the calculation, or we use up to 64 GPUs simultaneously to compute an ionisation yield curve.

Having performed some calculations with small and large basis sets for the same parameters, we cannot find any significant difference between them except the computation time. Thus all the calculations that we present are performed with a small basis set 6-31G.

3 Results and Discussions

This section is divided in two subsections. In the first one, we present some tests of the adequacy of our model in the case of atomic hydrogen exposed to a XUV pulse. Results for the water molecule are presented and discussed in the second subsection.

3.1 Verification of the model

In the case of atomic hydrogen, the atomic state wave functions are known analytically. It is therefore interesting to analyze to what extent the same wave function can be accurately reproduced in a basis of gaussian type orbitals (GTO) with the GAMESS(US) program. In Fig. 1, we compare the exact 1s-state wave function

(radial part) in momentum space with the corresponding wave function (radial part) generated by the GAMESS(US) program in two different bases of GTO’s, a small one (6-31G) and a bigger one (6-311++G**). A correlation consistent basis set d-aug-cc-pVTZ has also been used but the resulting radial part of the wave function which coincides with the one obtained with the 6-311++G** basis set, is not shown in the figure. These results show clearly that the GTO basis sets are not accurate in the region of small momenta. This is expected since small momenta correspond, in the configuration space, to large distances which are not well described by gaussians that are very localized functions. Although the small distances can be described relatively well using a small basis set (see the Figure 1), we have to account for the long range behaviour as well. In this case, a bigger basis set works better.

In Fig. 2a, we show the ionisation yield as a function of the frequency in the case where atomic hydrogen is exposed to a two-cycle pulse of zero carrier envelope phase and for 1014 W/cm2 peak intensity. The calculations have been performed in two ways: (i) by solving the TDSE with the Coulomb potential fully taken into account and (ii) by using our model in which we only take into account the ground state of atomic hydrogen. In the latter case, describing the atomic ground state with the exact wave function or in a basis of GTO’s does not affect the results despite the fact that GTO’s are unable to reproduce accurately the behaviour of the wave function for small momenta. We clearly see in Fig. 2a that the results obtained by solving TDSE and by using our model are very close to each other, although with a constant shift. To unveil the nature of this shift we performed two different lowest order perturbation theory (LOPT) calculations. In the first LOPT calculation, we considered single photon ionisation from the ground state of hydrogen atom to a continuum p-wave, described by a Coulomb wave. This is the LOPT TDSE curve, that coincides perfectly with the result of TDSE calculation. The second LOPT calculation has a plane wave as the final state, and everything else is the same as in the first case. This is the LOPT SPAM curve, that coincides perfectly with the SPAM calculation. We can conclude that the constant shift that we observe is due to the wrong description of the continuum, which is close to a plane wave in the case of SPAM model, while it is a Coulomb wave in reality. However, what is hidden behind the log-log scale of the plot is the fact that the absolute difference between the SPAM and TDSE becomes smaller and smaller. This shows that for high frequencies we can approximate a Coulomb wave with a plane wave. It turns out that using plane wave instead of Coulomb wave is equivalent to introducing a constant factor. We use this factor to correct the predictions for atomic hydrogen and for water molecule as well. In fact, this factor is not "ad hoc" as it can be derived from the dipole matrix elements in the case of a Coulomb wave and a plane wave final state of the ionising system. The comparison of the corrected SPAM result with full TDSE result is given in Figure 2b.

3.2 Water molecule

Let us now consider the water molecule. In this case, full TDSE calculations are not tractable and must rely on some approximation schemes. In the high frequency regime considered here and given the good agreement we obtained between TDSE and our model calculations for atomic hydrogen, it is reasonable to think that applying our approach to the water molecule considered as a single active electron system should give relatively good results.

In this subsection, we study the probability of orientation-dependent single ionization of the HOMO 1b1 by a laser pulse for various frequencies, peak intensities and carrier envelope phases. The reason why we chose water instead of a simpler multielectron molecule like H2 is that the water molecule is a very convenient multielectron quantum system for our purposes as: (a) its HOMO orbital is essentially a 2p atomic orbital of the oxygen atom with very little influence of each hydrogen atom (the value of magnetic quantum number for this 2p orbital depends on the particular position of the molecule in the reference frame); (b) ionisation of the water HOMO orbital leaves the geometry of the molecule unchanged, thus allowing us to apply the fixed nuclei approximation as the vibrational excitation is low.

Since the full final momentum wavepacket is available in the SPAM model, any information about the system can be easily extracted. The absence of the intermediate states means however that regimes where these states are important, like low frequency ionisation, cannot be treated accurately. Being a SAE approach one neglects all the dynamic interactions between the particles. Finally the Born Oppenheimer approximation neglects all the processes related to motion of the nuclei. Nevertheless, the SPAM model allows one to make predictions for any complex system, where the aforementioned approximations are adequate, in the single photon regime. The model is very scalable, so the limits on the size of the system are given by the hardware resources.

Fig. 3 shows two possible orientations of the water molecule with respect to the reference frame that are used later for the calculations. The origin of this reference frame is in the molecular plane. Let us recall that for all the cases treated here the electric field is polarized along the z-axis.

All the calculations have been performed in the 6-31G basis set. We are not aware of any significant discrepancies in our calculations that could be attributed to the incompleteness of this basis.

In Fig. 4, we consider the photon energy dependence of the probability of ionisation of the water molecule placed as in Fig. 3a, by a laser pulse with a fixed duration of 900 as, for three different angles between the polarization vector and the z-axis and for two peak intensities: $10^{15}$ W/cm2 and $5\times 10^{15}$ W/cm2. It can be seen that for a fixed peak intensity, the ionisation yield is the highest when the field polarization axis is along the axis of $2p_{z}$ oxygen orbital which practically coincides with the H2O HOMO orbital. The yield decreases rapidly when the polarization vector is rotated towards the x-axis. This conclusion supports the findings of Petretti et al.13. They showed that for a 800 nm wave length, the shape of the water molecule ionisation yield as a function of the orientation of the polarization vector is more or less the same as the shape of the molecular orbital from which the ionisation occurs. To stress this point, we show in Fig. 5, the ionisation yield

as a function of the angle $\theta$ between the polarization direction in the xz-plane and the z-axis. When the polarization axis is perpendicular to the z-axis (therefore in the xy-plane), we obtain an ionization yield which is equal to zero within the machine accuracy. This contrasts with Petretti’s results who obtained a small but finite yield for a frequency 10 times lower. Note that the ionisation yield gets flatter around $\theta=\pi/2$ for the lowest frequency, and also the overall shape shows significant variation.

In order to have some idea about how accurate is the prediction of the SPAM model for high frequencies, we ran the SAE-TDSE code from 13 in the 1-photon ionisation regime and compared the ionisation yield prediction of these two SAE models (see Figure 6). A significant difference between these models is the fact that for SAE-TDSE the full coulomb potential has been taken into account in generating the orbital basis, thus we have to correct the SPAM result with the same factor as for atomic hydrogen. Both ionisation yields are normalised to 1. In fact, we have perfect agreement between these models for high frequencies, and poor agreement for the photon energies near the ionisation threshold. This can be explained by the fact that SAE-TDSE uses 6000 Kohn-Sham orbitals to propagate the wavefunction, while in SPAM we don’t have any intermediate state at all. This discrepancy can be reduced as we include the lowest unoccupied orbital in the SPAM calculation, which will be done in future publications.

We can see in Figure 7 that the corrected SPAM coincides with SAE-TDSE in a wide intensity range as well. This agreement of the corrected SPAM model and SAE-TDSE approach in a wide intensity and frequency range indicates that the SPAM correction factor does not depend on intensity and on frequency.

4 Conclusions and perspectives

We have extended a new computationally inexpensive model, previously developed for atomic hydrogen to the treatment, within the SAE approximation, of the interaction of a complex quantum system with a high frequency ultrashort laser pulse. As a first application, we have applied this model to the single ionisation of the HOMO of the water molecule by an ultrashort XUV pulse. We studied the dependence of the single ionisation yield on the pulse frequency, the peak intensity, the orientation of the polarization vector and the carrier envelope phase. Our results clearly show that the model allows one to perform such calculations quickly for a complex system. Although it has been shown that this approach works rather well for atomic hydrogen in the low frequency regime we do not expect this model to be able to treat the ionisation of H2O by a low frequency laser pulse because the intermediate states for a many electron quantum system are very important. However we will try to investigate the behaviour of the model water molecule in a laser pulse having HOMO and LUMO orbitals included.

{acknowledgement}

The authors thank Prof. Clément Lauzin for many useful discussions. A.G. is "aspirant au Fonds de la Recherche Scientifique (F.R.S.-FNRS)". Yu.P. thanks the Université catholique de Louvain (UCL) for financially supporting several stays at the Institute of Condensed Matter and Nanosciences of the UCL. F.M.F. and P.F.O’M. gratefully acknowledge the European network COST (Cooperation in Science and Technology) through the Action CM1204 "XUV/X-ray light and fast ions for ultrafast chemistry" (XLIC) for financing several short term scientific missions at UCL. P.D. and A.G. acknowledge COST XLIC and F.R.S-FNRS for financing two short term scientific missions (STSM) in Trieste, Italy, and participation in COST XLIC meetings. The present research benefited from computational resources made available on the Tier-1 supercomputer of the Federation Wallonie-Bruxelles funded by the Region Wallonne under the Grant No. 1117545 as well as on the supercomputer Lomonosov from Moscow State University and on the supercomputing facilities of the UCL and the Consortium des Equipements de Calcul Intensif (CECI) en Federation Wallonie-Bruxelles funded by the F.R.S.-FNRS under the convention 2.5020.11. Y.P. is grateful to the Russian Foundation for Basic Research (RFBR) for financial support under the grant No. 16-02-00049-a.