High-order harmonic generation in atomic and molecular systems

TL;DR

This paper presents an analytical extension of the strong-field approximation to describe high-order harmonic generation in atomic and molecular systems, enabling direct study of target and laser pulse effects with good accuracy.

Contribution

It introduces an analytical method using a non-local separable potential to calculate key matrix elements in HHG, improving understanding of physical mechanisms and interference effects.

Findings

Analytical expressions match numerical solutions for atomic systems.

Method captures interference features in multicenter molecular systems.

Enables direct analysis of target and pulse influence on HHG spectra.

Abstract

High-order harmonic generation (HHG) results from strong-field laser matter interaction and it is one of the main processes that are used to extract electron structural and dynamical information about the atomic or molecular targets with sub-femtosecond temporal resolution. Moreover, it is the workhorse for the generation of attosecond pulses. Here we develop an analytical description of HHG, which extends the well established theoretical strong-field approximation (SFA). Our approach involves two innovative aspects: i) First, using a model non-local, but separable potential, we calculate the bound-free dipole and the rescattering transition matrix elements analytically for both atomic and molecular multicenter systems. In comparison with the standard approaches to the HHG process, these analytic derivations of the different matrix elements allows us to study directly how the HHG spectra depend on the atomic target and laser pulse features. We can turn on and of contributions having distinct physical origins or corresponding to different physical mechanisms. Our SFA results are compared, when possible, with the direct numerical integration of the time-dependent Schr\"odinger equation (TDSE) in reduced and full dimensionality. Excellent agreement is found for single and multielectronic atomic systems, modeled under the single active electron approximation, and for simple diatomic molecular systems. Our model captures also the interference features, ubiquitously present in every strong-field phenomenon involving a multicenter target.

Figures

Equations

$$\vec{\mu }(t)=-⟨\Psi (t)|\text{r}|\Psi (t)⟩,$$

$$\vec{\mu }(t)$$

$$I_{x{\rm N}}(\omega)\propto\>\Bigg{|}\int_{-\infty}^{\infty}{\textit{d}t\>e^{i\omega t}\>\vec{\mu}_{x{\rm N}}(t)}\Bigg{|}^{2}.$$

$$b(\text{p},t)=i{\int }_0^t\text{d}{\text{t}}^{\prime }\text{E}(t^{\prime })\cdot \text{d}\left[\text{p}+\text{A}(t^{\prime })\right]e^{-\text{i}S(\mathbf{p},t,t^{\prime })},$$

$$\textbf{d}_{\rm 1N}(\textbf{p}_{0})=\textit{i}\>\textbf{p}_{0}\>\frac{(p_{0}^{2}+\Gamma^{2})+(\frac{p_{0}^{2}}{2}+I_{p})}{(p_{0}^{2}+\Gamma^{2})^{\frac{3}{2}}(\frac{p_{0}^{2}}{2}+I_{p})^{2}}\Bigg{[}\frac{\Gamma+\sqrt{2I_{p}}}{2\pi\>(2I_{p})^{-1/4}}\Bigg{]}.$$

$${\vec{\mu }}_{1N}(t)=i{\int }_0^t\text{d}{\text{t}}^{\prime }\int {\text{d}}^3\text{p}$$

$${\vec{\mu }}_{1N}(t)$$

$$|\Psi (t)⟩$$

$${\vec{\mu }}_{2N}(t)$$

$$b(\text{p},t)$$

$${\mathbf{d}}_j({\mathbf{p}}_0)=-2\text{i}\mathcal{M}\frac{-{\text{p}}_0(3p_0^2+2I_p+2{\Gamma }^2)}{(p_0^2+{\Gamma }^2{)}^{\frac{3}{2}}(p_0^2+2I_p{)}^2}{e}^{-i{\text{R}}_j\cdot {\text{p}}_0}.$$

$${\vec{\mu }}_{2N}(t)=i\sum _{j=1}^2\sum _{j^{\prime }=1}^2$$

$${\vec{\mu }}_{2N}(t)$$

$${\vec{\mu }}_{11}(t)$$

$${\vec{\mu }}_{j{j}^{\prime }}(t)$$

$${\vec{\mu }}_{2N}(t)$$

$$I_{2N}(\omega )=I_{2N-\mathrm{Local}}(\omega )+I_{2N-\mathrm{Cross}}(\omega ),$$

$${\vec{\mu }}_{3N}(t)$$

$$b(\text{p},t)$$

$${\vec{\mu }}_{3N}(t)$$

$${\vec{\mu }}_{j{j}^{\prime }}(t)$$

$${\vec{\mu }}_{11}(t)$$

$${\vec{\mu }}_{3N-\mathrm{Local}}(t)={\vec{\mu }}_{11}(t)+{\vec{\mu }}_{22}(t)+{\vec{\mu }}_{33}(t),$$

$${\vec{\mu }}_{3N-\mathrm{Cross}}(t)={\vec{\mu }}_{3N-\mathrm{Cros}{\mathrm{s}}_1}(t)+{\vec{\mu }}_{3N-\mathrm{Cros}{\mathrm{s}}_2}(t),$$

$${\vec{\mu }}_{3N-\mathrm{Cros}{\mathrm{s}}_1}(t)={\vec{\mu }}_{12}(t)+{\vec{\mu }}_{21}(t)+{\vec{\mu }}_{23}(t)+{\vec{\mu }}_{32}(t),$$

$${\vec{\mu }}_{3N-\mathrm{Cros}{\mathrm{s}}_2}(t)={\vec{\mu }}_{13}(t)+{\vec{\mu }}_{31}(t).$$

$$\text{E}(t)={\mathcal{E}}_{0}f(t)\sin ({\omega }_{0}t+{\varphi }_0){\mathbf{e}}_z.$$

$${\omega }_{\text{cutoff}}=I_p+3.17U_p.$$

$$I_{2N-{\mathrm{R}}_1}(\omega )\propto {\left|{\int }_{-\infty }^{\infty }\text{d}t{e}^{i\omega t}[{\vec{\mu }}_{11}(t)+{\vec{\mu }}_{21}(t)]\right|}^2$$

$$I_{2N-{\mathrm{R}}_2}(\omega )\propto {\left|{\int }_{-\infty }^{\infty }\text{d}t{e}^{i\omega t}[{\vec{\mu }}_{22}(t)+{\vec{\mu }}_{12}(t)]\right|}^2.$$

Full text

I Introduction

High-order harmonic generation (HHG) is a conversion process resulting from the extremely high nonlinear interaction of a short and intense laser pulse with gas atoms or molecules or, recently, solid targets and nanostructures L’Huillier et al. (1991); Lyngå et al. (1996); Krausz and Ivanov (2009); Ghimire et al. (2011); Vampa et al. (2014); Ciappina et al. (in press). Nowadays, the HHG process is the conventional route for the production of spatially and temporally coherent extreme-ultraviolet (XUV) light, as well as light pulses in the sub-femtosecond and attosecond regimes Corkum and Krausz (2007). Coherent light sources in the ultraviolet (UV) to XUV spectral range are ubiquitously employed in a broad range of subjects, including basic research, material science, biology, and lithography Krausz and Ivanov (2009). Furthermore, the molecular HHG process encodes electronic orbital structure information and presents, as a consequence, a reliable method to retrieve molecular intrinsic parameters with attosecond and sub-Ångström temporal and spatial resolution, respectively Lein (2007); Lein et al. (2002a); Kanai et al. (2005); Torres et al. (2007); Morishita et al. (2008). Taking this objective in mind, several theoretical and experimental work have been conducted in order to optimize, improve and understand the molecular HHG process. Furthermore, HHG in atoms is one of the most studied topics of strong field physics and several theoretical models, besides of the solution of the time-dependent Schrödinger equation (TDSE), have been developed to describe it. Amongst them the most widely used and successful is the strong field approximation (SFA) Lewenstein et al. (1995, 1994).

The underlying physics of the HHG process is usually understood invoking the so-called “three step model”: (i) tunnel ionization; (ii) propagation in the laser field “continuum”, and (iii) recombination with the parent ion Corkum (1993); Lewenstein et al. (1994). According to this approach, when a strong laser pulse interacts with an atomic or molecular target a bounded electron is liberated through tunnel ionization (this happens when the laser electric field is close to its peak during an optical cycle). This “free” electron is then driven away from the ionic core and accelerated by the laser electric field, developing an oscillating trajectory. During this journey, the electron accumulates kinetic energy, that is released during the recombination process in the form of a high energy photon. As this three-step process usually occurs every half-cycle of the laser field, the spectrum of the generated coherent radiation consists of peaks at odd integer multiples of the driven laser frequency.

On the other hand, for multicenter molecules much less experimental Lyngå et al. (1996); Cornaggia et al. (1994); Schmidt et al. (1994); Lappas and Marangos (2000); Velotta et al. (2001); Hay et al. (2002); Lein et al. (2002b) and theoretical Ivanov and Corkum (1993); Zuo et al. (1993); Yu and Bandrauk (1995); Moreno et al. (1997); Kopold et al. (1998); Lappas and Marangos (2000); Alon et al. (1998); Averbukh et al. (2001); Bandrauk and Yu (1999); Kreibich et al. (2001) work have been done. The direct numerical solution of the TDSE for more complex systems with more than two centers is a quite challenging and formidable task from the numerical and computational viewpoints. Even for the simplest diatomic molecule, one has to solve a three-dimensional TDSE, that typically requires the utilization of a multicore CPUs and large amount of memory. In addition, the interpretation of the results extracted from TDSE is a not trivial task, in particular if one wants to disentangle the underlying physical mechanisms contributing to the total HHG spectrum.

As was mentioned above, the initial interest in the molecular HHG was due to the fact that it offers additional degrees of freedom and promising possibilities, such as the alignment of the molecular axis with respect to the laser field polarization axis. Specifically for a diatomic system, the existence of a distinctive quantum-interference minima pattern in the spectra and its dependence with the molecular orientation have been theoretically predicted by Lein et al. Lein et al. (2002a, 2003); Ciappina et al. (2007a). It is demonstrated that this pattern is due to a destructive interference from the high-harmonic emission at spatially separated centers and the internuclear distance can be accurately obtained scrutinizing the HHG spectra. In addition, the chance of controlling the phases and improve the phase-matching condition opens a route to the investigation in this area. More importantly, research on this field revealed how the distinctive features of the molecular HHG spectra can be used to retrieve structural information in simple molecules Itatani et al. (2004). Furthermore, the high harmonic generation spectroscopy has shown the possibility to extract structural and dynamical information from the molecular HHG spectra in more complex targets (for a couple of examples see e.g. Ciappina et al. (2007b, 2008); Laarmann et al. (2007)). Finally, studies in small molecules shows that the temporal evolution of the electronic wavefunction can be directly recovered Smirnova et al. (2009); Haessler et al. (2010); Kraus et al. (2015).

In this paper, we use the SFA within the framework of the Lewenstein’s model to study the HHG from atomic, diatomic and three atomic molecular systems in the few-cycle IR laser pulse regime. The derivation for the two and three centers molecular systems is constructed as a consecutive extension of the atomic model. For simplicity, our analytical model is based on a non-local potential which is approximately a short-range (SR) potential. We compute HHG spectra for those three systems and, for the atomic case, compare the results with the numerical solution of the 3D-TDSE in the single active electron (SAE) approximation.

Our approach involves two innovative aspects: i) First, using a model non-local, but separable potential, we calculate the bound-free dipole and the rescattering transition matrix elements analytically for both atomic and molecular multicenter systems. In comparison with the standard approaches to the HHG process, these analytic derivation of the different matrix elements allows us to study directly how the HHG spectra depend on the atomic target and laser-pulse features; we can turn on and off contributions having distinct physical origins or corresponding to different physical mechanisms. This allow us to quantify their weights in the various regions of the HHG spectra; ii) Second, as in Ref. Suárez et al. (2016), in our theory the dipole matrix elements are free from non-physical gauge and coordinate system dependent terms – this is achieved by adapting the coordinate system, in which SFA is performed, to the centre from which the corresponding part of the time dependent wave function originates. We compare, when possible, our SFA results with the numerical solutions of the time-dependent Schrödinger equation in full dimensionality. Excellent agreement is found for atomic and molecular systems, including multielectronic systems modeled under the SAE. Our model captures also the interference features, ubiquitously present in any multicenter target.

This article is organized as follows. In Section II, we address the main theory that describes the HHG process and, in particular, the derivation of the time-dependent dipole matrix element within the SFA. Here, we make use of the results previously presented in Suárez et al. (2015, 2016) to obtain the analytical expressions needed to compute each of the individual contributions. Particularly for the three-center molecular system, we develop a new set of equations to compute the time-dependent dipole matrix element, making use of the non-local short-range potential bound states. In section III, HHG spectra for the atomic and molecular cases are numerically calculated. Results in hydrogen and argon atoms are presented, comparing them with those obtained from the 3D-TDSE. For diatomics, we analyze two systems: H${}_{2}^{+}$ and H2. The basic analysis of the interference minima of the harmonic spectra with respect to the alignment for H${}_{2}^{+}$ is discussed. In addition, the contribution of the different processes to the total spectra is assessed. A time analysis of the HHG spectra using a Gabor transformation is performed and the influence of the short and long trajectories is investigated. In addition, CO2 and H2O define our three-center molecular systems. For both cases, we investigate the dependence of the HHG spectra with respect to the molecular orientation and extract information about the different mechanism contributing to the total HHG spectra. Finally, in section IV, we summarize the main ideas and present our conclusions. Atomic units will be uses throughout the manuscript otherwise stated.

II Theory of HHG within the SFA

In this section we develop a quantum mechanical approach of HHG using the generalized SFA model described in Ref. Lewenstein et al. (1995, 1994). The source of the additional frequencies that are generated during the interaction of a strong laser pulse with an atomic or molecular target, is the nonlinear dipole oscillation of the medium. Therefore, the aim is to calculate this dipole response by mean of the solution of the time dependent Schrödinger equation. The time-dependent radiation dipole moment reads:

[TABLE]

where $|\Psi(t)\rangle$, is the state describing the time-evolution of the atomic or molecular system under study. In general, within the SFA statement, we can write the wavefunction of the whole system as a superposition of the ground, $|0\rangle$, and continuum states, $|\textbf{v}\rangle$, as $|\Psi(t)\rangle=e^{\textit{i}I_{p}\textit{t}}(a(t)|0\rangle+\>\int{\textit{d}^{3}\textbf{v}\>\textit{b}(\textbf{v},t)|\textbf{v}\rangle})$, where the transition amplitude of the continuum states is denoted by $\textit{b}(\textbf{v},t)$. After some algebra with the above equations and only considering transitions from the bound to the continuum states the time-dependent dipole radiation moment reads

[TABLE]

where the bound-continuum transition dipole matrix element is defined as $\textbf{d}(\textbf{v})=-\langle\textbf{v}|\textbf{r}|0\rangle$. The radiation emitted by a single atom is proportional to the time-dependent dipole moment $\vec{\mu}(t)$. In this way the harmonic spectrum, $I(\omega)$, is calculated as the modulus squared of the Fourier transformed dipole acceleration -$\textbf{a}(t)$- related to the defined time-dependent dipole matrix element, Eq. (2), by the Ehrenfest theorem as $|\tilde{\textbf{a}}(\omega)|=|\omega^{2}\>\tilde{\mu}(\omega)|$. In this way we can compute the harmonic spectra as:

[TABLE]

In here the subscript $x$ represents the total numbers of atoms in the system to study and it will count as $x=1,2,\dots,n$, where $n$ is the total number of atoms of the molecule. For case of a diatomic molecule, i.e. constituted by two atoms, the subscript reads as: $x=2$ meaning $\vec{\mu}_{\rm 2N}(t)$ and $I_{\rm 2N}(\omega)$. Notice that both the atomic and molecular harmonic spectra depend directly on the time-dependent dipole moment which in turn depends on the form of the bound-continuum matrix element and the continuum states transition amplitude, that are different for each of atomic, diatomic or multiatomic system under study.

II.1 Calculation of the time-dependent dipole moment for atomic systems: $\vec{\mu}_{\rm 1N}(t)$

In order to have all the ingredients to compute the harmonic spectrum for an atomic system, $I_{1N}(\omega)$, using the Eqs. (2) and (3) we need to know the exact dependency of the bound-continuum matrix element and the continuum states transition amplitude. The method to find the transition amplitude of the continuum states and bound-continuum matrix element for an atom under the influence of an intense laser pulse has been described in our previous work Suárez et al. (2015). We therefore take advantage of those results and only explain here the new derivations needed to tackle the HHG problem.

The transition amplitude for the continuum states of the atomic system reads as:

[TABLE]

where the exponent phase factor is “the semiclassical action” ${S}({\bf p},t,t^{\prime})=\int_{t^{\prime}}^{t}{\>d{\tilde{t}}\left\{[{\bf p}+\textbf{A}({\tilde{t}})]^{2}/2+I_{p}\right\}}$ defining all the possible electron trajectories from the birth time $t^{\prime}$ until the “recombination” one $t$.

The explicit expression for the bound-continuum transition dipole matrix element obtained in Suárez et al. (2015) is

[TABLE]

Inserting Eqs. (4) and (5) in the time-dependent dipole moment, Eq. (2), and changing variables to the canonical momentum defined by $\textbf{p}=\textbf{v}-\textbf{A}(t)$ we get,

[TABLE]

Equation (6) has to be understood as follows: the electron is ionized at time $t^{\prime}$ with a certain probability defined by, $\textbf{E}(t^{{}^{\prime}})\cdot\textbf{d}_{\rm 1N}\left[\textbf{p}+\textbf{A}(t^{{}^{\prime}})\right]$. During its excursion in the continuum the electronic wavefunction is then propagated until the time $t$ acquiring a classical phase $\textit{S}(\textbf{p},t,t^{\prime})$ to finally recombine with the ion core at time $t$ with a rate given by $\textbf{d}_{\rm 1N}^{*}[\textbf{p}+\textbf{A}(t)]$. All possible combinations of birth time and momenta must be considered and therefore a multidimensional integration is required, where their contributions are added up coherently. Note that Eq. (6) configures a highly oscillatory integral, both in the momentum p and $t^{\prime}$ variables. As a consequence it is convenient to rewrite the integral over ${\bf p}$ using the stationary-phase approximation or saddle point method. In order to do that is necessary to find the extremal points over the exponential phase. The extrema $\mathbf{p}={\bf p}_{s}$ are found from the solutions of $\nabla_{{\bf p}}\textit{S}(\textbf{p})|_{{\bf p}_{s}}={\bf 0}$. These saddle point momenta ${\bf p}_{s}$ thus can be written as $\textbf{p}_{s}=-\frac{1}{\tau}\int_{t^{{\prime}}}^{t}{\textbf{A}(\tilde{t})\textbf{d}\tilde{t}}$. Here, $\tau=t-t^{\prime}$ is the excursion time of the electron in the continuum. Expanding, the function $\textit{S}(\textbf{p},t,t^{\prime})$ in a Taylor series around the roots $\textbf{p}_{s}$ and then applying the standard saddle point method to the momentum integral over ${\bf p}$, the time-dependent dipole matrix element for the atomic system reads as:

[TABLE]

where we have introduced an infinitesimal parameter, $\varepsilon$, to avoid the divergence at $t=t^{\prime}$ (for a detailed discussion see Suárez et al. (2015, 2016)). The harmonic spectrum, $I_{1N}(\omega)$, is then numerically computed inserting Eq. (7) in Eq. (3).

II.2 Calculation of the time-dependent dipole moment for diatomic molecular systems: $\vec{\mu}_{\rm 2N}(t)$

In order to calculate the harmonic spectrum generated by a diatomic molecule we use the results obtained in Ref. Suárez et al. (2016). As we can extract from that reference the general wavefunction describing the state for a diatomic molecule can be written as:

[TABLE]

from which the molecular time-dependent dipole moment $\vec{\mu}_{\rm 2N}(t)$ is easily obtained and have the following form:

[TABLE]

In the above equation we require to insert the explicit expression for the continuum states transition amplitude $b(\textbf{p},t)$:

[TABLE]

and the bound-continuum dipole matrix element ${\bf d}_{j}$.

In the derivation of the length-gauge SFA model for HHG in diatomic molecules, and in particular for the computation of the bound-continuum dipole matrix element $\textbf{d}(\textbf{v})=-\langle\textbf{v}|\textbf{r}|0\rangle$, an unphysical term is neglected, without give a consistent reason/argument (see Hetzheim et al. (2007); Chirilă and Lein (2006); Busuladžić et al. (2008) for more details). This term, that is a linear function of the internuclear distance R, immediately introduces convergence problems as ${\textbf{R}}\to\infty$. Clearly, this behavior introduces conflicts between the length and velocity gauges predictions, observed in the case of above-threshold ionization (ATI) as well. The root of the problem relies in the degree of approximation to handle the continuum states, considered as a set of plane waves for the molecular system, without considering the relative position of each atomic core. This creates an unphysical treatment and therefore the appearance of such unphysical term. In our approach we solve this issue by computing ${\bf d}_{j}(\textbf{v})=-\langle\textbf{v}|(\hat{\textbf{r}}-\textbf{R}_{j})|0_{j}\rangle$, where here the bound-continuum dipole matrix element is calculated with respect to each atomic center located at $\textbf{R}_{j}$. Note that if no approximations are done, i.e. if we consider the case where $\langle\textbf{v}|$ is a scattering wave of the field free Hamiltonian H0, the above mentioned problem will not arise – the scattering waves are orthonormal to the ground states $|0_{j}\rangle$. However, as the main core of the SFA is to handle the continuum states as Volkov states, i.e. neglecting the influence of the residual molecular potential once the electron is in the continuum, the convergence problems would remain if we do not correct the bound-continuum dipole matrix element. The full derivation of the bound-continuum dipole matrix element for the ATI problem in a two-center molecular system was introduced in Ref. Suárez et al. (2016) - this bound-continuum dipole matrix element is the same used for the computation of HHG. In addition, an extended derivation for a three-center molecular system is presented in the Appendix A.

The bound-continuum dipole matrix element is defined by $\textbf{d}_{\rm 2N}(\textbf{p}_{0})=\sum_{j}^{2}{\bf d}_{j}({\bf p}_{0})$, where ${\bf d}_{j}$ denotes the bound-continuum dipole matrix element related to the nucleus located at the position ${\bf R}_{j}$ and is given by:

[TABLE]

Here $\mathcal{M}$ is a normalization constant (for details see Ref. Suárez et al. (2016)). Note that the index $j$ can take the value of 1 (or 2), referring to the nucleus located at the position ${\bf R}_{1}$ (or ${\bf R}_{2}$) on the left (and right).

The time-dependent radiation dipole moment $\vec{\mu}_{\rm 2N}(t)$ thus reads:

[TABLE]

where subscript $j$ and $j^{\prime}$ represent the ionization the recombination atom positions, respectively.

Equation (12) contains information about all the recombination processes occurring in the entire molecule during the HHG phenomenon and can then be written as a sum of components as:

[TABLE]

The four terms in the above equation encode all the individual molecular recombination processes. Our physical interpretation of those contributions is as follows:

1. (i)

An electron is ionized from the atom placed at the Left with respect to the coordinate origin at time $t^{\prime}$ with certain probability: $\textbf{E}(t^{\prime})\cdot\textbf{d}_{1}\left[\textbf{p}+\textbf{A}(t^{\prime})\right]$. During its excursion in the continuum this electron accumulates a phase which depends on the position from where it was detached, in this case $\textbf{R}_{1}$. Finally, because the electric field changes its sign and the electron returns to the parent ion, the probability of recombination results $\textbf{d}^{*}_{1}[\textbf{p}+\textbf{A}(t)]$. In this step the excess of energy acquired from the laser electric field is converted into a high energy photon. This whole process is described by:

[TABLE] 2. (ii)

The second term is understood in a similar way. In this case the ionization and recombination processes occur in the core placed at the Right. The equation describing this process, $\vec{\mu}_{22}$, is similar to Eq. (14) but considering the dipole matrix element $\textbf{d}_{2}$ and now the electron is detached from the position $\textbf{R}_{2}$. The two processes described before are spatially localized (involving only one core placed at a fixed position $\textbf{R}_{1}$ or $\textbf{R}_{2}$) and we then refer to them as “Local Processes”. 3. (iii)

The last two terms, $\vec{\mu}_{21}(t)$ and $\vec{\mu}_{12}(t)$, describe events involving two atoms at two different positions $\textbf{R}_{1}$ or $\textbf{R}_{2}$. Here $\vec{\mu}_{21}$ can be understood as follows: the electron is tunnel-ionized from the atom on the Right with certain probability given by: $\textbf{E}(t^{\prime})\cdot\textbf{d}_{2}\left[\textbf{p}+\textbf{A}(t^{\prime})\right]$. After this ionization event the electron starts to move under the laser electric field influence accumulating energy and acquire a phase: $e^{-i\left\{S({\bf p},t,t^{\prime})+\textbf{R}_{2}\cdot\left[\textbf{A}(t)-\textbf{A}(t^{\prime})\right]\right\}}$. Finally the electron returns back to the other core (Left) at the time $t$ to end up its journey in a recombination process that has an amplitude proportional to: $\textbf{d}^{*}_{1}[\textbf{p}+\textbf{A}(t)]$. As in previous cases the excess energy is emitted in a form of a high energy photon. Considering both centers are involved in the HHG process, we call these terms as “Cross processes”. The equation describing these processes reads:

[TABLE]

where now $j\neq j^{\prime}$ denotes the nucleus-index located at left ($j$=1) or right ($j$=2).

Note from the above description that we have to account four different possible processes corresponding to four different time-dependent transition dipole moments. Two of them are “Localized” and two “ Cross” representing all the possible recombination scenarios in our diatomic molecule.

Similarly to the atomic case, in order to obtain the molecular time-dependent dipole matrix $\vec{\mu}_{\rm 2N}(t)$ we apply the saddle point method in the momentum variable p. In fact, the phases of the local contributions in Eq. (14), function on the relative positions $\textbf{R}_{1/2}$ of the atoms, cancel each other defining a saddle point momentum ${\bf p}_{s}$ equivalent to the one presented in the atomic case (see Sec. II.A). On the other hand, the cross process presents more complex phases, that directly depend on the position variables. For instance, in the process $\vec{\mu}_{21}\>(\>\vec{\mu}_{12})$ the saddle-point momentum can be found to be: $\textbf{p}_{s+}=-\frac{1}{\tau}\left[\textbf{R}+\int_{t^{\prime}}^{t}{\textbf{A}(\tilde{t})\textbf{d}\tilde{t}}\right]\left(\textbf{p}_{s-}=-\frac{1}{\tau}\left[-\textbf{R}+\int_{t^{\prime}}^{t}{\textbf{A}(\tilde{t})\textbf{d}\tilde{t}}\right]\right)$. In all our cases we are going to work with shorter internuclear distances, where the following condition is fulfilled: $R<\mathcal{E}_{0}/\omega^{2}$, with $\mathcal{E}_{0}$ and $\omega_{0}$ being the laser electric field peak amplitude and carrier frequency, respectively. As a consequence, it is not needed to consider this saddle-point momentum definition (for more details about the validity of this approximation see Chirilă and Lein (2006)). Thus, we proceed by applying the standard saddle point momentum to all the local and cross contributions. The total time-dependent dipole moment for our diatomic molecule then reads as:

[TABLE]

Finally the total HHG spectrum can be calculated using Eq. (3), similarly to the atomic case, but using the time-dependent dipole matrix obtained in Eq. (16). As it was discussed, four terms are needed to compute each molecular harmonic spectrum. Each term represents a different process and this is equivalent to the split made in the time-dependent dipole matrix element. We label each contribution depending on the position of the atoms, e.g. from the $\textit{Left}-\textit{Left}$ term we obtain the $I_{\rm 2N,11}(\omega)$ spectrum. Similarly we write the other three terms as $I_{\rm 2N,22}(\omega)$, $I_{\rm 2N,12}(\omega)$ and $I_{\rm 2N,21}(\omega)$, respectively.

It is convenient to identify two main contributions in the total harmonic spectrum, Eq. (16), namely, (i) one generated for the Local processes and (ii) other developed by the Cross processes. In this way we can write the total harmonic spectrum as:

[TABLE]

where $I_{\rm 2N-Local}(\omega)=I_{\rm 2N,11}(\omega)+I_{\rm 2N,22}(\omega)$ and $I_{\rm 2N-Cross}(\omega)=I_{\rm 2N,12}(\omega)+I_{\rm 2N,21}(\omega)$ denote the local and cross terms, respectively.

II.3 Calculation of the time-dependent dipole moment for three-center molecular systems: $\vec{\mu}_{\rm 3N}(t)$

The computation of the HHG spectrum generated by a three-center molecule using the definition in Eq. (3) involves the search of the exact bound states describing the whole system. In order to do so we use a method similar to the one presented in Ref. Suárez et al. (2015, 2016). In short, we consider a three-center molecule as a set of three atoms placed at $\textbf{R}_{1}=-\frac{\textbf{R}}{2}$, $\textbf{R}_{2}=0$ and $\textbf{R}_{3}=\frac{\textbf{R}}{2}$, respectively, where R is the so-called internuclear distance, defined as the distance between the atoms placed at $\textbf{R}_{1}$ and $\textbf{R}_{3}$, when the molecule is linear. The state describing the time evolution of a three-center molecule can be written as a coherent superposition of the states $|\Psi(t)\rangle$ as $|\Psi(t)\rangle=e^{\textit{i}I_{p}\textit{t}}\left(a(t)|0\rangle+\sum_{j=1}^{3}\int{\textit{d}^{3}\textbf{v}\>\textit{b}_{j}(\textbf{v},t)|\textbf{v}\rangle}\right)$, where the subscript $j=1,2,3$, refers to the contributions of the spatially localized nuclei at $\textbf{R}_{1}$, $\textbf{R}_{2}$ and $\textbf{R}_{3}$, respectively. By employing the Schrödinger equation on that state and our basic SFA approach, the molecular time-dependent dipole moment $\vec{\mu}_{3N}(t)$ reads:

[TABLE]

$\vec{\mu}_{\rm 3N}(t)$ is defined as a superposition of the bound-continuum dipole matrix of each atom on the molecule, i.e. $\textbf{d}_{3N}(\textbf{v})=\sum_{j=1}^{3}\textbf{d}_{j}(\textbf{v})$. The exact dependency of the bound-continuum matrix element is presented in the Appendix A (see Eq. (63) for more details).

Using the exact definition of the bound-continuum matrix element, the total continuum states transition amplitude, $\textit{b}(\textbf{v},t)=\sum_{j=1}^{3}b_{0,\mathit{j}}(\textbf{p},t)$, reads as:

[TABLE]

The explicit expression for the molecular time-dependent dipole matrix element $\vec{\mu}_{3N}(t)$ is obtained inserting Eq. (19) in Eq. (18). As in the case of diatomics, it is also possible here to disentangle each of the recombination processes contributing to the total harmonic spectrum. In order to do so we write $\vec{\mu}_{\rm 3N}(t)$ as a sum of nine terms as:

[TABLE]

The above equation contains information about all the possible recombination scenarios present in our three-center molecule. In order to make clearer the interpretation let us write the individual time-dependent dipole matrix element $\vec{\mu}_{jj^{\prime}}(t)$ explicitly as

[TABLE]

where the subscripts “$j$” and “$j^{\prime}$” refer to the position $\textbf{R}_{1}$, $\textbf{R}_{2}$ and $\textbf{R}_{3}$ of each of the atoms in the three-center molecule. In Eq. (21) the first subscript $j$ represents the atom from where the electron is detached and can be $j=1,2,3$. In addition, the second one, $j^{\prime}$, labels the atom where the recombination process occurs, and can also take the values 1,2 or 3.

Note that, as in the case of atoms and diatomic molecules, in Eq. (21) we have applied the saddle point method in the momentum p integral. In Eq. (21) we have followed the same criteria as in the diatomic system (see Sec. II.B) and in this way the saddle point momentum $\mathbf{p}_{s}$ is the conventional one.

As in the case of diatomics, the nine terms of Eq. (20) represent the Local and Cross processes. These different terms should be understood as follows:

1. (i)

The first term, $\vec{\mu}_{11}$, describes the process of an electron ionized from the atom placed at $\textbf{R}_{1}$ at time $t^{\prime}$ with probability: $\textbf{E}(t^{\prime})\cdot\textbf{d}_{1}\left[\textbf{p}_{s}+\textbf{A}(t^{\prime})\right]$. This electron, during its excursion in the continuum, accumulates a phase which depends on the position from where it was detached, in this case $\textbf{R}_{1}$. Finally, because the change in the sign of the laser electric field, the electron returns to the parent ion, with a recombination probability given by $\textbf{d}^{*}_{1}[\textbf{p}_{s}+\textbf{A}(t)]$. As a result of this recombination stage the excess of energy acquired from the laser electric field is converted into a high energy photon. As an example the time dependent dipole equation describing this process, where $j=1$ and $j^{\prime}=1$, reads as:

[TABLE] 2. (ii)

The second and third terms, $\vec{\mu}_{22}$ and $\vec{\mu}_{33}$, describe the same process, but for the atoms located at $\textbf{R}_{2}$ and $\textbf{R}_{3}$, respectively. These three process are spatially localized: the electron starts and ends at the same point, the same ion core. We then refer to them as “Local processes”. 3. (iii)

From the fourth to the seventh terms we have the cross processes with the closer neighbor in one and other direction. In this case notice that in our reference frame the second atom is placed at $\textbf{R}_{2}=0$. These processes are understood as in the diatomic cases. 4. (iv)

The last two terms are also cross processes. For instance, the eighth term can be understood as follows: one electron tunnels ionize from the atom located at $\textbf{R}_{1}$ with probability: $\textbf{E}(t^{\prime})\cdot\textbf{d}_{1}\left[\textbf{p}_{s}+\textbf{A}(t^{\prime})\right]$. This electron starts to move with the electric field and acquires a phase $e^{-i\left\{S({\bf p}_{s},t,t^{\prime})+\textbf{R}_{1}\cdot\left[\textbf{A}(t)-\textbf{A}(t^{\prime})\right]\right\}}$. It then recombines at time $t$ with the farthest ion-core at $\textbf{R}_{3}$ with an amplitude $\textbf{d}^{*}_{3}[\textbf{p}_{s}+\textbf{A}(t)]$. The last term is understood in a similar way, but inverting the tunnel ionization and recombination positions.

For our three-center molecular system is also possible to group the processes as Local and Cross. As in the diatomic case the sum of all these terms represents the total time-dependent dipole element, $\vec{\mu}_{3N}(t)=\vec{\mu}_{\mathrm{3N-Local}}(t)+\vec{\mu}_{\mathrm{3N-Cross}}(t)$. In the same way we can split the contributions depending on the excursion of the electron in the continuum before recombination. The shorter excursions are represented by the local processes where only one atom is involved. For the cross processes we have two possibilities: the recombination with (i) the closest neighbor or (ii) with the farthest one. Those contributions are denoted by:

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

In order to describe the direct processes we set $j=j^{\prime}$. For instance, the Local process for the Right atom located at $\textbf{R}_{1}$ is described by the time-dependent dipole matrix element $\vec{\mu}_{11}(t)$. On the other hand, the Cross processes are those where $j\not=j^{\prime}$.

Finally, in order to compute the total time-dependent dipole element $\vec{\mu}_{\rm 3N}(t)$ of our three-center molecule we need to evaluate each of the contributions defined by Eq. (21). The HHG spectrum can then be obtained by Fourier transforming $\vec{\mu}_{\rm 3N}(t)$ (see Eq. (3)). The separation of the time-dependent dipole matrix element allows us to compute the harmonic spectrum from each process separately as we will see in the next Section.

III Results and Discussion

In this section we calculate the harmonic spectra for different systems using the equations previously presented. In addition, we compare the harmonic spectra emitted from hydrogen and argon atoms with the exact numerical solution of the 3D-TDSE. A scan over different laser wavelength and peak intensities is performed in order to verify and validate the model. In a second stage, we apply our molecular approach to two prototypical diatomic systems: H${}_{2}^{+}$ and H2. We display the different contributions coming from the local and cross recombination processes which helps to distinguish which contributions are interfering constructively and destructively to the total high harmonic spectra. Finally, we present results for more complex molecules: CO2 and H2O. For these cases, besides to disentangle the different contributions to the HHG spectra, we analyze the influence of the angular orientation.

The numerical integration of Eqs. (7), (16) and (21) has been performed by employing a rectangular rule with dedicated emphasis on the convergence of the results. The HHG process is driven by an ultrashort laser pulse with an electric field of the form:

[TABLE]

The field has a carrier frequency $\omega_{0}=\frac{2\pi c}{\lambda_{0}}$, where $c$ is the speed of light ($c\approx 137$ a.u) and $\lambda_{0}$ the central laser wavelnegth and $\mathcal{E}_{0}$ is the field peak amplitude, linearly polarized in the $z$-axis. ${f(t)}=\sin^{2}(\omega_{0}t/2N_{c})$ denotes the pulse envelope, with $N_{c}$ the total number of cycles, and the parameter $\phi_{0}$ is the carrier envelope phase (CEP). Under the dipole approximation, the influence of the magnetic field component of the laser field is neglected.

III.1 Atomic systems. Comparison between SFA and 3D-TDSE models

To calculate the harmonic spectra of an hydrogen atom we perform the Fourier transform of the time-dependent dipole moment presented in Eq. (7). We set $\Gamma=1$ and $\gamma=38$ a.u in our non-local potential in order to match the ionization potential $I_{p}=0.5$ a.u. of the hydrogen atom. We consider a pulse with $N_{c}=4$ total cycles and $\phi_{0}=0$ rad. A total of $131072$ points in the time window $t\in[0,t_{\rm F}]$, where $t_{\rm F}=N_{c}T_{0}$ and $T_{0}=2\pi/\omega_{0}$, are used during the numerical integration. The simulation of the harmonic spectra for H at different laser wavelengths and using our quasiclassical SFA model is shown in Fig. 1(a). In addition, in Fig. 1(b) we show the HHG spectra obtained by using the numerical solution of the 3D-TDSE.

In order to compute the HHG spectra displayed in Fig. 1 we consider the laser pulse described by Eq. (27), with a laser peak intensity of $I_{0}=1.58\times 10^{14}$ W$\,\cdot$ cm*-2* and different laser wavelengths (see the panels label for details). In order to calculate the HHG spectra of Fig. 1(b) we numerically solve the 3D-TDSE in the length gauge. Thus, by Fourier transform of the dipole acceleration, calculated from the time propagated electronic wave function, the HHG spectra is obtained. We have used our code, which is based on an expansion of spherical harmonics, $Y_{lm}$ considering only the $m=0$ terms due to the cylindrical symmetry of the problem. The numerical technique is based on a Crank-Nicolson method implemented on a splitting of the time-evolution operator that preserves the norm of the wave function.

Both panels of Fig. 1 reveal the typical HHG behavior, namely (i) a rapidly decreasing of the harmonic yield for the lower harmonic orders ($<10^{\mathrm{th}}$); (ii) a plateau with almost constant yield and (iii) an abrupt end at the so-called HHG cutoff. The cutoff energy is one of the most important features of any HHG spectrum. It can be defined as the maximum photon energy that can be released at recollision. Classically it is possible to prove that Corkum (1993); Lewenstein et al. (1994):

[TABLE]

where $\omega_{\textrm{cutoff}}$ is the maximum photon energy and $U_{p}=I_{0}/4\omega_{0}^{2}$ is the ponderomotive energy. As can we see from Fig. 1 both the SFA and 3D-TDSE calculations show the expected classical cutoff defined by Eq. (28), noted with a dashed line of each color at $\omega_{\textrm{cutoff}-800}=1.59$ a.u. (43.26 eV), $\omega_{\textrm{cutoff}-1200}=2.97$ a.u. (80.8 eV) and $\omega_{\textrm{cutoff}-1400}=3.87$ a.u. (105.3 eV), respectively. From Eq. (28) we should note that $\omega_{\textrm{cutoff}}\propto I\lambda^{2}$ and this behaviour can also be observed in Fig. 1. For instance, the spectra at $\lambda_{3}=1400$ nm have a cutoff energy about 4 times higher than the one calculated using a wavelength of $\lambda_{1}=800$ nm.

A natural next step would be to test our model with a more complex atom. In order to do so in Fig. 2 we show HHG spectra for an argon atom, calculated both with by (i) our quasiclassical SFA (Fig. 2(a)) and (ii) using the numerical solution of the 3D-TDSE under the SAE approximation (Fig. 2(b)). We employ two different intensities and using a laser pulse with a central frequency of $\omega_{0}=0.057$ a.u., that corresponds to a wavelength of about 800 nm. As in the previous case, we confirm that our model is capable to capture not only the dependency of the harmonic spectra with the wavelength, but also with the laser peak intensity. As we can see, and considering that $I_{2}>I_{1}$, a clear cutoff extension in the HHG spectra for $I_{2}$ is observed. A remarkable good agreement between both methods is clearly seen in Fig. 2 and for both laser intensities.

The HHG spectra presented both for a single electron system (H, Fig. 1) and a complex target (Ar, Fig. 2) reveal the excellent agreement between our quasiclassical SFA model and the numerical solution of the 3D-TDSE.

III.2 Diatomic molecular systems

In this section we calculate HHG spectra for two prototypical diatomic molecules: H${}_{2}^{+}$ and H2.

III.2.1 H${}_{2}^{+}$ molecule

Figure 3 shows the numerically computed HHG spectra for an H${}_{2}^{+}$ molecule by using the quasiclassical SFA model presented in Sec. II.B. The H${}_{2}^{+}$ system is modeled by two identical centers separated by an internuclear distance $R=2.2$ a.u. (1.16 Å) and the molecular axis forms a $\theta$ angle with respect to the incident laser electric field polarization, i.e. $\textbf{R}=(0,0,R\cos\theta)$. The parameters of our non-local potential are set to $\Gamma=1.0$ and $\gamma=0.1$ a.u. in order to reproduce the minimum at the equilibrium internuclear distance, $R_{0}$ = 2.0 a.u. (1.06 Å), in the potential energy surface (PES). In our short-range potential toy model the total ionization potential extracted from the potential energy surface yields $I_{p}=0.68$ a.u. (18.50 eV). We compute this electronic ground state energy to fix the asymptotic behavior of the H${}_{2}^{+}$ potential energy surface (see Ref. Suárez et al. (2016) for more details). This last value differs from the one obtained with a real Coulomb potential that leads a pure electronic energy of $1.1$ a.u. (30 eV) approximately.

The incident laser field shape is identical to the one used in the atomic case and has a central frequency $\omega_{0}=0.057$ a.u., corresponding to a wavelength $\lambda=800$ nm and photon energy of $1.55$ eV, respectively. The total number of cycles is $N_{c}=4$ -this defines a full-width at half-maximum FWHM value of $5.2$ fs- and $\phi_{0}=0$ rad. The time step is set to $\delta t=0.032$ a.u. and this corresponds to a total of $N_{t}=20000$ points for the numerical integration. The time window is $t\in[0,t_{\rm F}]$, where $t_{\rm F}\approx 11$ fs denotes the final time, i.e. the end of the laser electric field pulse. Finally, the laser peak intensity is set to $I_{0}=5\times 10^{14}$ W$\,\cdot$ cm*-2*.

In Fig. 3 we display results for a scan of four different molecular orientations, namely Fig. 3(a) $\theta=0^{\circ}$ (this value corresponds to the so-called parallel alignment), Fig. 3(b) $20^{\circ}$, Fig. 3(c) $40^{\circ}$ and Fig. 3(d) $45^{\circ}$, respectively. As we can see in all the panels an absolute minimum over the total HHG spectra is clearly visible and the harmonic order where these minima are located increases with the orientation angle. The existence of those minima and their dependency with the alignment angle can be explained by invoking an interference phenomenon as we will see below. In the most simplest picture the minima appears as a consequence of the harmonic emission of two radiant points (see Lein et al. (2002a) for more details).

According to the equation describing the destructive interference of two radiant sources: $R\cos\theta=(2m+1)\lambda_{k}/2$, where $\lambda_{k}$ is the de Broglie wavelength of the returning electron and considering the “fundamental” instance $m=0$, the minima should be located at the $18^{\mathrm{th}}$, $20^{\mathrm{th}}$, $30^{\mathrm{th}}$ and $36^{\mathrm{th}}$ harmonic order for $\theta=0^{\circ}$, $\theta=20^{\circ}$, $\theta=40^{\circ}$ and $\theta=45^{\circ}$, respectively. The positions of the minima for our SFA calculation are $\approx 35^{\mathrm{th}}$, $\approx 37^{\mathrm{th}}$, $\approx 45^{\mathrm{th}}$ and $\approx 54^{\mathrm{th}}$, respectively (see the vertical lines in all the panels of Fig. 3). We speculate that the shifts in harmonic frequency are related with the kind of potential used in our calculations; the short range potential does not correctly describe the low energy part of the HHG spectra, where the Coulomb potential plays an important role Ciappina et al. (2007a). We note, however, that our SFA calculation for $\theta=40^{\circ}$ is in excellent agreement with the numerical solution of the 2D-TDSE and 3D-TDSE for the H${}_{2}^{+}$ molecule Lein et al. (2002a, 2003). Lastly, we observe that in all the HHG spectra of Fig. 3 the position of the classical cutoff is in excellent agreement with Eq. (28) (see the arrows in all the panels of Fig. 3). Particularly, for our H${}_{2}^{+}$ molecular system and the laser parameters used in our simulations, the cutoff frequency is $\omega_{\textrm{cutoff}}=4.15$ a.u. (112.92 eV), corresponding to the $72^{\mathrm{th}}$ harmonic order.

Clearly, our quasiclassical molecular SFA model has drawbacks and advantages. The first advantage is from the computational viewpoint; the numerical calculations using the our SFA model are much faster than the numerical solution of the 3D-TDSE. The computation of one single HHG spectrum for a set of fixed parameters takes few seconds. The second, and might be the most important one, is the possibility to disentangle the different components contributing to the final harmonic spectra (see Sec. II.B). In order to do so in the Fig. 4 we plot the different contributions to the HHG spectra for an H${}_{2}^{+}$ molecule aligned at $\theta=20^{\circ}$ with respect to the laser field polarization. Figure 4(a) shows the main contributions of the harmonic spectra, namely the total $I_{2N}(\omega)$ (red circle line) , the local $I_{\rm 2N-Local}(\omega)$ (blue line) and the cross, $I_{\rm 2N-Cross}(\omega)$ (dark brown asterisk line) (for details see Sec. II.B). As we can see from this picture the two-center destructive interference is not present in neither in the local nor in the cross contributions. The latter have a deep minimum but at a different position, about the $60^{\mathrm{th}}$ harmonic order, while the local contribution remains almost constant in yield for all the harmonic frequencies. In order to trace out the origin of the two-center destructive interference present in the total HHG spectra in Fig. 4(b) we plot the contributions depending on the recombination atom, calculated as:

[TABLE]

and

[TABLE]

From this figure we can clearly see that there is a deep minimum for both terms and it is located at the same position as for the total HHG spectra. It means that, for the case of the recombination on $\textbf{R}_{1}$ (dark green circle line), the electron wavepacket ionized at $\textbf{R}_{1}$ interferes with the one coming from $\textbf{R}_{2}$ and the other way around. These minima are then generated by the destructive interference of such electron wavepackets. From the drawbacks side, we have seen that our short range non-local potential is unable to accurately reproduce the interference minima positions for some of the molecular orientation angles. We note, however, that these minima are typically washed out when an average over the molecular orientation is considered, configuration that is commonly used in molecular HHG experiments.

III.2.2 Time-frequency analysis for H${}^{+}_{2}$

We have seen in Fig. 4 that the independent processes $\vec{\mu}_{11}(t)/\vec{\mu}_{22}(t)$ and $\vec{\mu}_{21}(t)/\vec{\mu}_{12}(t)$, are the ones interfering and creating the deep minimum in the total HHG spectra. In order to dig deeper about the existence of this distinctive feature a Gabor analysis Chirilă et al. (2010); Gabor (1946) over the different contributions is displayed in Fig. 5. The Gabor transformation was performed upon the time-dependent dipole moment calculated using our quasiclassical SFA model. The laser parameters are the same as in Fig. 4.

This time-frequency analysis allows us to reveal the half-cycle bursts of radiation from which the HHG spectrum is composed and the main trajectories contributing. In Fig. 5(a) and Fig. 5(b) we show the local and cross process, $\vec{\mu}_{11}(t)$ and $\vec{\mu}_{21}(t)$, respectively. As we can observe from these figures, they both look almost equal and similar to the atomic case. In both cases we have the contribution of the short and long trajectories. For the earlier cycles, the 1st and 2nd, the short trajectory contributions dominate while for the latest cycles both long and short trajectories have the same weight. The main differences between these two contributions are in the low energy region around the end of the laser pulse, the 3rd optical cycle, where the contribution of the cross processes $\vec{\mu}_{21}(t)$ is slightly higher than the local ones.

Finally, we plot the mixed, Fig. 5(c), and total, Fig. 5(d), contributions. From these figures is evident the presence of an interference minimum for the whole temporal window. This means that the, $\vec{\mu}_{11}(t)$ and $\vec{\mu}_{21}(t)$ processes, that describe electrons arriving at the same point $\textbf{R}_{1}$ from two different atomic sources, $\textbf{R}_{1}$ and $\textbf{R}_{2}$ respectively, cancel each other and an interference zone is seen for an harmonic order of around $35^{\mathrm{th}}$. These two contributions are dominated by the short-trajectories, therefore both incoming electron wavepackets arrive at the same time and as a consequence a destructive interference is observed. This feature is inherited to the total time-dependent dipole moment (see Fig. 5(d)).

III.2.3 H2 molecule

The next simplest diatomic molecule is the H2. In order to investigate the behaviour and versatility of our semiclassical SFA model, we compute HHG spectra using the time-dependent dipole moment presented in Sec. II. We consider an H2 molecule in equilibrium where the two H atoms are separated a distance of $R=1.4$ a.u. (0.74 Å) The ionization potential of the outer electron predicted by our non-local potential is $I_{p}=1.5$ a.u. (40.82 eV) and it was calculated by setting $\Gamma=1.0$ and $\gamma=0.12$ a.u. With these parameters our model reproduce the PES of H2 with a minimum at the equilibrium internuclear distance Lein (2007). The driven laser pulse has the same parameters as for the case of H${}_{2}^{+}$.

Figure 6 shows the different contributions to the total HHG spectrum by considering two different molecular orientations: parallel, $\theta=0^{\circ}$, (Fig. 6(a)) and perpendicular, $\theta=90^{\circ}$, (Fig. 6(b)) with respect to the incident laser pulse polarization. The total HHG spectra (in red) is computed as the sum of all possible processes (see Sec. II.B for details). In both panels we have grouped two main contributions: (i) the Local and (ii) the Cross ones. The Local contributions (blue line) are processes related with only one atom or position, i.e. they involve the sum of processes involving only one single atom, meaning ionization from the $\textbf{R}_{1}/\textbf{R}_{2}$ and recombination in the same atom. On the other hand, the Cross contributions (in dark brown) include processes involving both of the atoms in the molecule, i.e ionization from the atom located at $\textbf{R}_{1}$ and recombination on the atom located at $\textbf{R}_{2}$ and the other way around.

The first observation regarding Fig. 6 is that for the case of parallel orientation, Fig. 6(a), the total HHG spectrum starts to gradually decrease for harmonic orders higher than the $\approx 30^{\mathrm{th}}$. This behaviour is due to a destructive interference of the local and the cross processes. The latter shows a deep minimum around the $\approx 40^{\mathrm{th}}$ order. In the case of the molecule perpendicularly oriented, Fig. 6(b), an extended plateau with a cutoff around $90^{\mathrm{th}}$ is clearly visible. In both cases, parallel and perpendicular, the molecular HHG spectra shows a deep minimum around the $12^{\mathrm{th}}$ harmonic order. As in the case of H${}_{2}^{+}$ previously presented, the utilization of a short range potential restricts our results to the higher order harmonics, where the influence of the molecular potential details is less relevant.

It is interesting to note that for the case of perpendicular orientation, $\theta=90^{\circ}$, (Fig. 6(b)), both the local and cross terms contribute evenly in the plateau region of the HHG spectra, while for the parallel orientation $\theta=0^{\circ}$, (Fig. 6(a)) the local and cross contributions present a different behavior. We can then infer that for the $\theta=90^{\circ}$ case the total harmonic spectrum reaches a maximum yield at the cutoff region. This is due to the fact that, for this favourable orientation, the contribution of each of the processes, local and cross, is comparable.

Finally, in Fig. 7, we show the total HHG spectra for three different molecular orientations, $\theta=0^{\circ}$, $\theta=45^{\circ}$ and $\theta=90^{\circ}$ and an averaged case over nine values of $\theta$ in the range $[0^{\circ}-360^{\circ}]$. Our diatomic molecule is symmetrical with respect to the origin, i.e. $\textbf{R}_{1}=-\textbf{R}/2$ and $\textbf{R}_{2}=\textbf{R}/2$ and, consequently, the total HHG spectra for $\theta=0^{\circ}$ and $\theta=180^{\circ}$ are identical. The same behaviour is observed for the spectra at $45^{\circ}$, $135^{\circ}$, $225^{\circ}$ and $315^{\circ}$ or for $90^{\circ}$ and $270^{\circ}$. We can observe in Fig. 7 how different molecular configurations contribute to the total HHG spectra. As we can see the difference in the total spectra for different orientation angles is hardly to notice for lower harmonic orders ($<30^{\mathrm{th}}$). Differences start to be noticeable in the mid-plateau and cutoff regions. In these spectra ranges the larger HHG yield is reached for the perpendicular orientation ($\theta=90^{\circ}$), thus confirming the results presented in Fig. 6. Two final remarks are in order, namely (i) the averaged total harmonic spectra is about one order of magnitude lower than the one at perpendicular orientation; (ii) the average procedure washes out any two-slit interference fingerprint.

III.2.4 Time-frequency analysis for H2

In Fig. 8 we perform a Gabor transformation upon the time-dependent dipole moment for both an H atom and our diatomic H2 molecule. Our aim with this time-energy analysis is to investigate the influence of the short and long trajectories for the molecular system and highlight the differences with the atomic case. In Fig. 8(a) we show the calculation for the H atom, while in Fig. 8(b) we show the same analysis for the molecular system randomly oriented. In both cases we have considered a laser peak intensity of $I_{0}=5\times 10^{14}$ W$\,\cdot$ cm*-2*.

In general both figures look quite different. The atomic system, Fig. 8(a), is mostly dominated by the short trajectories while the molecular system, Fig. 8(b), have a prevailing contribution from the long ones. This is so because the orientation average procedure removes any fingerprint of the molecular interferences.

From a detailed comparison between the atomic and molecular cases we observe that for the former, even when the short trajectories are dominant at the beginning of the laser pulse (first 2 optical cycles), some contribution of the long ones survives for the later optical cycles where long and short contributions are similar (3 optical cycle). On the contrary, in the molecular system short and long trajectories contribute to different optical cycles. For instance, in the 1st and 2nd optical cycles the main contribution is from the short trajectories while for the 3nd and 4rd optical cycles a big contribution of the long trajectories appears. In the molecular system the contributions of the long trajectories start to increase, being paramount for the later optical cycles where the contribution of the short one is less significant.

III.3 Three-center molecular systems

In order to study systems with more degrees of freedom and describe the different processes contributing to the total HHG spectra, as we have done for diatomics, we apply the model described in Sec. II.C to both CO2 and H2O molecules.

III.3.1 The carbon dioxide molecule

The carbon dioxide molecule CO2 is a linear system formed by three atoms, O$=$C$=$O, where the two oxygen atoms are separated by a distance of $R=4.38$ a.u (2.31 Å) when the system is in equilibrium. At this equilibrium state the parameters of the non-local potential are set to $\Gamma=0.8$ and $\gamma=0.11$ a.u corresponding to an ionization potential of $I_{p}=0.50$ a.u. (13.6 eV). This value is in excellent agreement with the actual CO2 ionization energy (13.77 eV) CO2 (2016). The incident laser electric field is defined in Eq. (27) and we use a laser wavelength and peak intensity of $\lambda=800$ nm and $I_{0}=1\times 10^{14}$ W$\,\cdot$ cm*-2*, respectively. The laser pulse has four total cycles (11 fs of total duration) and the CEP is set to $\phi_{0}=0^{\circ}$.

The HHG spectra computed by using our quasiclassical SFA model for the CO2 system is presented in Fig. 9. In the Fig. 9(a) we show the different contributions to the HHG spectra: the total $I_{\rm 3N}(\omega)$ (solid line with red circles), calculated from the time-dependent dipole presented in Eq. (20), the local $I_{\rm 3N-Local}(\omega)$ (blue solid line), computed with Eq. (23), and the cross $I_{\rm 3N-Cross}(\omega)$ (dark brown line with asterisks), extracted from Eq. (24). These calculations show the well known HHG plateau that ends with a cutoff (marked with a red arrow) at around the $21^{\mathrm{th}}$ harmonic order (this last value is in perfect agreement with the one predicted by the semiclassical law Eq. (28)). Both local and cross contributions have almost the same yield over all the frequency range and only minor differences are visible.

In Fig. 9(b) we present a split of the local processes, $I_{3N-11}(\omega)$ (solid line with purple circles), $I_{3N-22}(\omega)$ (solid line with light blue squares) and $I_{3N-33}(\omega)$ (dashed line). As we can see the contribution from the O atoms, placed at the end of the molecule, is equal in amplitude and shape and different in yield from the $I_{3N-22}(\omega)$ (corresponding to the C atom placed at the origin). This means that the O atoms contribute a slightly less than the C atom placed at the origin. We notice, however, that the shapes and positions of the minima are the same for the three contributions.

In Fig. 9(c) we present each of the contributions that build up the total cross processes. We have separated them depending on how long the ionized electron travels in the continuum before recombination. The Cross1 (solid line with orange circles), Eq. (25), and the Cross2 (solid brown line), Eq. (26) have similar yields. The main difference between these two HHG spectra is the yield: the cross1 has a bigger contribution than the Cross2. The position of the absolute minima around the $19^{\mathrm{th}}$ harmonic order is present in both contributions as in the local term. For the calculations in Figs. 9(a)-(c) we consider a CO2 molecule aligned perpendicular to the incident laser pulse polarization, i.e. the internuclear axis vector is forming a angle of $\theta=90^{\circ}$, being this the most favorable configuration (see Fig. 9(d)).

Finally, in Fig. 9(d) we present a set of total HHG spectra for different molecular orientations, namely parallel ($\theta=0^{\circ}$), oblique ($\theta=45^{\circ}$) and perpendicular ($\theta=90^{\circ}$). In addition we include an averaged HHG spectrum, obtained coherently adding # orientations. We can observe a similar behavior as for the case of H2 (see Fig. 7), i.e. the difference in the total spectra for different orientation angles is hardly to see for lower harmonic orders and starts to be visible in the mid-plateau and cutoff regions. Furthermore, the perpendicular orientation appears to be the dominant one. The comparable behavior between the CO2 and H2 molecules support the fact that the former could be considered as a ’stretched’ diatomic O2 molecule for interference minima calculations Lein (2007).

III.3.2 The water molecule

One of the most important three-center molecules is water (H2O) since it is part of the building blocks of biological life. In this section we theoretically investigate harmonic spectra of the H2O molecule using our semiclassical SFA approach.

We consider an H2O molecule under the influence of the strong laser field described by Eq. (27). The H2O molecule is an angular molecule with two H atoms and one O atom. At equilibrium the internuclear distance of the bond H$=$O is about $R=1.8$ a.u. (0.95 Å) and the angle between the two H atoms $\alpha=104.5^{\circ}$. For these parameters, and considering an ionization potential of $I_{p}=0.46$ a.u. (12.52 eV)H2O (2016), we set the parameters of our non-local potential to $\Gamma=0.8$ and $\gamma=0.1$ a.u.

In Fig. 10 we show HHG spectra for a laser wavelength and peak intensity of $\lambda=800$ nm and $I_{0}=1\times 10^{14}$ W$\,\cdot$ cm*-2*, respectively. The laser pulse has four total cycles (11 fs of total duration) and the CEP is set to $\phi_{0}=0$ rad. In Fig. 10(a) we show HHG spectra for both five different molecular orientations, $\theta=0^{\circ}$, $\theta=20^{\circ}$, $\theta=45^{\circ}$, $\theta=60^{\circ}$ and $\theta=90^{\circ}$ and an averaged case. The molecular axis is fixed in space and forms an angle of $\alpha/2$ with respect to the vector position $\textbf{R}_{1}$. Furthermore, $\theta$ defines the angle between this molecular axis and the laser electric field polarization (see Fig. 12 and the Appendix A for more details).

The dependency of the HHG spectra with respect to the molecular orientation is quite evident. For lower harmonic orders all the orientations appear to be equivalent and the main differences start to materialize for the harmonic orders $\gtrsim 12^{\mathrm{th}}$. As we can see in Fig. 10(a) the HHG spectra for $\theta=0^{\circ}$ (solid line with asterisks), $20^{\circ}$ (solid line with left-pointing triangle) and $45^{\circ}$ (dashed line) exhibit a similar structure. The other two orientations, $60^{\circ}$ (right-pointing triangle line) and $\theta=90^{\circ}$ (square line), present an harmonic yield several orders of magnitude lower in this region. The total HHG spectra for all the molecular orientations show a slight minimum around the $17^{\mathrm{th}}$ harmonic order that could be attributed to interference effects, although it is not an easy task to characterize it using a simple interference formula as in the case of diatomics.

We have also included in Fig. 10(a) an averaged HHG spectra over 8 values of $\theta$ in the range $[0^{\circ}-360^{\circ}]$ (dashed red line). As we can see the minimum survives the orientation average. Furthermore, for $\theta=90^{\circ}$ (square line) the total HHG spectra rapidly decreases for harmonic orders $>16^{\mathrm{th}}$. This means that the interference between the local and cross processes is destructive and function of the molecular orientation. This behavior introduces a decrease of the total HHG yield. We note that for H2O, contrarily to the CO2 case, an enhancement of the total HHG spectra is observed when the molecule is oriented parallel, $\theta=0^{\circ}$, to the laser electric field polarization. As we have done both for diatomics and CO2 in Fig. 10 (b) we plot the different terms contributing to the total HHG averaged spectra. Contrarily to the oriented case, here the local and cross processes appear to constructively contribute to the total HHG spectrum.

In order to study more deeply the underlying physics behind the enhancement and decrease of the total HHG spectra for $0^{\circ}$ and $90^{\circ}$ we plot in Fig. 11 the different contributions for these two cases. For the case of $\theta=0^{\circ}$, Fig. 11(a), the decrease of the HHG yield is evident for harmonic orders higher than the $15^{\mathrm{th}}$. Around this harmonic order both contributions, the local and cross, have a similar yield and the coherent sum develops in a destructive interference decreasing the total HHG spectra in about 3 orders of magnitude. On the other hand, for $\theta=90^{\circ}$, Fig. 11(b), we observe a steadily decrease of the cross processes, of about two order of magnitude, in the whole spectral range. Consequently, we can argue that in this case the cross contribution is almost negligible (solid brown line with squares) and the total HHG spectra is dominated by the local processes (solid blue line).

IV Conclusions and Outlook

We present a quasiclassical approach that deals with molecular HHG within the SAE. Our model could be considered as a natural extension to the one introduced for above-threshold ionization in atoms  Suárez et al. (2015) and molecules Suárez et al. (2016). The focus of our study is on di- and tri-atomic molecular systems, although the extension to more complex systems appears to be straightforward. Firstly, we have validated our formalism comparing the atomic HHG spectra with results extracted from the 3D-TDSE and using a large set of laser intensities and wavelengths. For the molecular systems we have shown our approach is able to capture the interference features, ubiquitously present in every molecular HHG. The main advantages of our model are: (i) the possibility to disentangle the underlying contributions to the HHG spectra. In this way, we could isolate local and cross processes and also treat both fixed and randomly oriented molecules; (ii) the low computational cost. By considering our approach involves only 1D and 2D time integrations, all the other quantities are analytical, it is clear that we compute molecular HHG spectra without too much computational effort; (iii) the concrete feasibility to model complex molecular ground states. For all the studied molecular cases we were able to model reasonable well the initial molecular ground state, varying the parameters of our non-local potential.

Acknowledgements.

This work was supported by the project ELI–Extreme Light Infrastructure–phase 2 (CZ.02.1.01/0.0/0.0/15_008/0000162 ) from European Regional Development Fund, Ministerio de Economía y Competitividad through Plan Nacional (FIS2011-30465-C02-01, FrOntiers of QUantum Sciences (FOQUS): Atoms, Molecules, Photons and Quantum Information FIS2013-46768-P, FIS2014-56774-R, and Severo Ochoa Excellence Grant SEV-2015-0522), the Catalan Agencia de Gestio d’Ajuts Universitaris i de Recerca (AGAUR) with SGR 2014-2016, Fundació Privada Cellex Barcelona and funding from the European Union’s Horizon 2020 research and innovation Programme under the Marie Sklodowska-Curie grant agreement No. 641272 and Laserlab-Europe (EU-H2020 654148). N.S. was supported by the Erasmus Mundus Doctorate Program Europhotonics (Grant No. 159224-1-2009-1-FR-ERA MUNDUS-EMJD). N. S., A. C. and M. L. acknowledge ERC AdG OSYRIS and EU FETPRO QUIC. J.B. acknowledges FIS2014-51478-ERC and the National Science Centre, Poland – Synfonia grant 2016/20/W/ST4/00314. J.A.P.-H. acknowledges support from Laserlab Europe (Grant No. EU FP7 284464) and the Spanish Ministerio de Economía y Competitividad (FURIAM Project No. 78 FIS2013-47741-R and PALMA project FIS2016-81056-R)

Appendix A Strong field approximation for three-center molecular systems

In this section we develop an analytical model to obtain the direct probability transition amplitude, as well as the bound and scattering states, necessary to calculate the HHG spectra for three-center molecular systems. This approach can be considered an extension of the atomic and diatomic models presented in Refs. Suárez et al. (2015, 2016). Our quasiclassical formalism takes advantage of both the single active electron and dipole approximations. In addition, we consider the nuclei of the molecule are fixed in space (the so-called frozen core approximation).

Direct transition probability amplitude

We consider a molecular system of three independent atoms, as is shown in Fig. 12, under the influence of an intense and short laser field. In the limit when the wavelength of the laser, $\lambda_{0}$, is larger compared with the Bohr radius, $a_{0}=0.529$ nm, the electric field of the laser beam around the interaction region can be considered as spatially homogeneous. This means that the interacting atoms do not experience any spatial dependence of this driving field. Then, only its time-variation is taken into account (the above asseverations define the so-called dipole approximation).

Therefore, the laser electric field can be written as:

[TABLE]

The field has a carrier frequency $\omega_{0}=\frac{2\pi c}{\lambda_{0}}$, where $c$ is the speed of light, $\mathcal{E}_{0}$ the laser electric field peak amplitude, and we consider the laser field is linearly polarized along the $z$-direction. In Eq. (31) ${f(t)}$ denotes the envelope of the laser pulse and $\phi_{0}$ defines the CEP (see Sec. III for more details).

The TDSE that describes the whole laser-molecule interaction can be written as:

[TABLE]

where $\hat{H}_{0}=\frac{{\hat{\bf p}}^{2}}{2}+\hat{V}(\textbf{r})$ defines the laser-field free Hamiltonian, with ${\hat{\bf p}}=-i{\nabla}$ the canonical momentum operator and $\hat{V}(\textbf{r})$ the potential operator that describes the interaction of the nuclei with the active electron. $\hat{V}_{int}({\bf r},t)=-q_{e}\hat{{\bf E}}(t)\cdot\hat{{\bf r}}$ represents the interaction of the molecular system with the laser radiation, written in the dipole approximation and length gauge. $q_{e}$ denotes the electron charge (in atomic units $q_{e}=-1.0$ a.u.).

We shall restrict our model to the low ionization regime, where the SFA is valid Keldysh (1965); Faisal (1973); Reiss (1980); Grochmalicki et al. (1986); Lewenstein et al. (1994, 1995). Therefore, we work in the tunneling regime, where the Keldysh parameter $\gamma=\sqrt{I_{p}/2U_{p}}$ ($I_{p}$ is the ionization potential of the system and $U_{p}=\frac{\mathcal{E}_{0}^{2}}{4\omega_{0}^{2}}$ the ponderomotive energy acquired by the electron during its incursion in the laser field) is less than one, i.e. $\gamma<1$. In addition, we assume that $V(\textbf{r})$ does not play an important role in the electron dynamics once the electron appears in the continuum.

These observations, and the following three statements, define the standard SFA, namely:

1. (i)

Only the ground, $|0\rangle$, and the continuum states, $|\textbf{v}\rangle$, are taken into account in the interaction process. 2. (ii)

There is no depletion of the ground state $(U_{p}<U_{sat})$. 3. (iii)

The continuum states are approximated by Volkov states; in the continuum the electron is considered as a free particle solely moving driven by the laser electric field.

For a more detailed discussion of the validity of the above statements see e.g. Refs. Lewenstein et al. (1994, 1995); Suárez et al. (2015).

Based on (i), we propose a state, $|\Psi(t)\rangle=\sum_{j=1}^{3}|\Psi_{j}(t)\rangle$, to describe the time-evolution of the three-center system, i.e. a superposition of three atomic states. In turn, each independent state, $|\Psi_{j}(t)\rangle$, is a coherent superposition of ground, $|0\rangle=\sum_{j=1}^{3}|0_{j}\rangle$, and continuum states, $|\textbf{v}\rangle$  Lewenstein et al. (1994, 1995):

[TABLE]

where the subscript $j=1,2,3$ refers to the position $\textbf{R}_{1}$, $\textbf{R}_{2}$ and $\textbf{R}_{3}$ of each of the atoms in the three-center molecule, respectively.

The factor, $a(t)$, represents the amplitude of the ground state which is considered constant in time, $a(t)\approx 1$, under the assumption that there is no depletion of the ground state. The latter follows directly from statement (ii). The pre-factor, $e^{\textit{i}I_{p}\textit{t}}$, represents the phase oscillations which describe the accumulated electron energy in the ground state ($I_{p}=-E_{0}$ is the ionization potential of the molecular target, with $E_{0}$ the ground state energy of the three-center molecular system). Furthermore, the transition amplitudes to the continuum states are denoted by $\textit{b}_{j}(\textbf{v},t)$, with $j=1,2,3$ depending on the atomic nuclei. These amplitudes depend both on the kinetic momentum of the outgoing electron and the laser pulse. Therefore, our main task is to derive a general expression for each transition amplitude ${b}_{j}({\bf v},t)$. In order to do so, we substitute Eq. (33) in Eq. (32). We shall consider that $\hat{H_{0}}|0_{1,2,3}\rangle=-I_{p}|0_{1,2,3}\rangle$ and $\hat{H_{0}}|\textbf{v}\rangle=\frac{\textbf{v}^{2}}{2}|\textbf{v}\rangle$ fulfill for the bound and continuum states, respectively. Consequently, the evolution of the transition amplitude ${b}_{j}({\bf v},t)$ becomes:

[TABLE]

Note that we have assumed that the electron-nucleus interaction is neglected once the electron appears at the continuum, i.e. $V({\bf r})|{\bf v}\rangle=0$, which corresponds to the statement (iii). Therefore, by multiplying Eq. (34) by $\langle{\bf v}^{\prime}|$ and after some algebra, the time variation of the transition amplitude $\textit{b}_{j}(\textbf{v},t)$ reads:

[TABLE]

The first term on the right-hand of Eq. (35) represents the phase evolution of the electron in the oscillating laser electric field. In the second term we have defined the bound-free transition dipole matrix element as:

[TABLE]

The state $|\textbf{v}\rangle$ represents a scattering state constructed as a plane wave, $|{\bf v}_{p}\rangle$ plus corrections on each center position, $|\delta{\bf v}_{j}\rangle$.

Based on statement (iii) our formulation only considers the continuum state as a plane wave $|{\bf v}_{p}\rangle$ for the calculation of the bound-free dipole matrix element. We shall pay special attention to the computation of Eq. (36). Notice that the plane waves are not orthogonal to the bound states due to the fact that the latter are defined depending on the relative position of each of the atoms, $\textbf{R}_{j}$, with respect to the origin of coordinates (see Suárez et al. (2016) for more details). So, for the $\textbf{R}_{j}$ contribution we introduce a correction to the dipole matrix element as:

[TABLE]

The third term of Eq. (35) defines the continuum-continuum transition matrix element. In our case we are interested in to describe processes where the electron is ionized and goes to the continuum, never returning to the vicinity of the remaining ion core, i.e. the so-called direct processes. As the direct ionization process should have a larger probability compared with the continuum-continuum one Lewenstein et al. (1995); Suárez et al. (2015), one might neglect the last term in Eq. (35), $\langle\textbf{v}|\textbf{r}|\textbf{v}^{\prime}\rangle\approx\textbf{0}$.

This is what we refer as zeroth order solution:

[TABLE]

The latter equation is easily solved by conventional integration methods (see e.g. Elsgoltz (2003)) and considering the Keldysh transformation Ehlotzky (1992); Keldysh (1965). Therefore, the solution can be written as:

[TABLE]

Notice that $j=1,2,3$ represents either an atom located at $\textbf{R}_{1}$, $\textbf{R}_{2}$ or $\textbf{R}_{3}$, respectively. For instance, to obtain the transition amplitude for the atom placed at $\textbf{R}_{1}$ we need to set $j={1}$ in Eq. (39).

Note that the above equation is written in terms of the canonical momentum ${\bf p}={\bf v}-{\bf A}(t)$  Lewenstein et al. (1994). Here, we have considered that the electron appears in the continuum with kinetic momentum ${\bf v}(t^{\prime})={\bf v}-{\bf A}(t)+{\bf A}(t^{\prime})$ at the time $t^{\prime}$, where v is the final kinetic momentum (note that in atomic units ${\bf p}={\bf v}$), and $\textbf{A}(t)=-\int^{t}{\textbf{E}(t^{\prime})dt^{\prime}}$ is the associated vector potential.

Equation (39) has a direct physical interpretation: it can be understood as the sum of all the ionization events that occur from the time $t^{\prime}$ to $t$. Then, the instantaneous transition probability amplitude of an electron at a time $t^{\prime}$, at which it appears into the continuum with momentum ${\bf v}(t^{\prime})=\textbf{p}+\textbf{A}(t^{\prime})$, is defined by the argument of the $[0,t]$ integral in Eq. (39). Furthermore, the exponent phase factor denotes the “semi-classical action”, ${S}_{j}({\bf p},t,t^{\prime})$, that defines a possible electron trajectory from the birth time $t^{\prime}$, at position $\textbf{R}_{j}$, until the “recombination” one $t$ as:

[TABLE]

Note that the transition amplitude equations obtained so far depend on the position from which the electron is tunnel ionized to the continuum. The semi-classical action ${S}_{j}({\bf p},t,t^{\prime})$ contains this dependency as well.

Considering we are interested in to obtain the transition amplitude $b_{0,j}({\bf p},t)$ at the end of the laser pulse, the time $t$ is set at $t=t_{\rm F}$. Consequently, we shall define the integration time window as $t\in[0,t_{\rm F}]$. Furthermore we set ${\bf E}(0)={\bf E}(t_{\rm F})={\bf 0}$, in such a way to make sure that the laser electric field is a time oscillating wave and does not contain static components (the same arguments apply to the vector potential ${\bf A}(t)$). Finally, the total transition amplitude for the direct process taking place on our three-center molecular system reads as:

[TABLE]

Bound states calculation

In this section, we are going to develop analytical expressions to obtain the bound states for our three-center molecular system. As in the atomic and diatomic cases we have chosen a non-local potential to describe the interaction of the electrons with the nuclei. We consider a molecular system with three fixed nuclei under the single active electron approximation. Our purpose is to find the analytical dependency of the bound state wavefunctions, that allow us to compute the bound-continuum matrix transition dipole and the direct transition amplitudes, Eq. (41).

The hamiltonian for the molecular system in the momentum representation can be written in a similar way as for the diatomic case, i.e. :

[TABLE]

where the first term on the right-hand side is the kinetic energy operator and the second one describes the interacting non-local potential between the active electron and each molecular nuclei:

[TABLE]

where $\phi({\bf p})=\frac{1}{\sqrt{{\bf p}^{2}+\Gamma^{2}}}$ is an auxiliary function and $\gamma^{\prime}=\frac{\gamma}{3}$ is a parameter related with the shape of the ground state. By using $\hat{H}_{\textbf{M}}(\textbf{p},\textbf{p}^{\prime})$ from Eq. (42), we write the stationary Schrödinger equation as follows:

[TABLE]

where $E_{0}$ denotes the energy of the bound state. Thus, for our three-center system the Schrödinger equation reads:

[TABLE]

Defining new variables $\check{\varphi}_{j}$ as:

[TABLE]

we could analytically obtain the bound states by solving Eq. (45) in the momentum representation. Explicitly we can write:

[TABLE]

where $I_{p}$ denotes the ionization potential that is related to the ground state potential energy by $E_{0}=-I_{p}$. The wavefunction $\Psi_{0\textbf{M}}(\textbf{p})$ for the bound state in momentum space is defined as:

[TABLE]

In order to find the value of the constants we multiply and divide Eq. (48) by $e^{-\textit{i}\textbf{R}_{j}\cdot\textbf{p}}$ and $\sqrt{p^{2}+\Gamma^{2}}$, respectively. After some algebra we find that:

[TABLE]

where $I_{1}$, $I_{2}$ and $I_{3}$ read as:

[TABLE]

[TABLE]

and

[TABLE]

respectively.

Solving the system of equations Eq. (49) we find the relations between the $\check{\varphi}_{j}$ defined by Eqs. (46) and $I_{1}$, $I_{2}$ and $I_{3}$. The system Eq. (49) is solved with the restriction:

[TABLE]

and

[TABLE]

$I_{1}$, $I_{2}$ and $I_{3}$ are written in spherical coordinates as:

[TABLE]

[TABLE]

and

[TABLE]

Finally, Eq. (48), can be written as:

[TABLE]

where $\mathcal{M}_{3N}=\gamma^{\prime}\>\check{\varphi}_{1}=\frac{\gamma}{3}\>\check{\varphi}_{1}$ is a normalization constant. It can be calculated using the usual normalization condition:

[TABLE]

From the above equation ${\mathcal{M}_{3N}}$ can be written as:

[TABLE]

with $I_{4}$ defined as:

[TABLE]

With the exact knowledge of $\mathcal{M}_{3N}$ we have now defined the bound state in our three-center molecular system from Eq. (58). Notice that the dependency of the system energy with the internuclear distance appears from the solution of the system of equations, Eq. (49). In order to find potential energy surface describing the whole molecule we need to use Eq. (53) to obtain the relation:

[TABLE]

Bound-continuum transition matrix element

The total dipole matrix element for the three-center molecular system is defined as a sum:

[TABLE]

Equation (63) can be explicitly written as:

[TABLE]

where:

[TABLE]

Finally notice that we could extract the contributions of each center from Eq. (64), i.e.

[TABLE]

[TABLE]

and

[TABLE]