Probing the role of excited states in ionization of acetylene

TL;DR

This study models acetylene ionization using VUV laser pulses and finds that excited states significantly influence ionization efficiency, harmonic generation, and photoelectron spectra, depending on laser wavelength and molecular alignment.

Contribution

It demonstrates the impact of excited states on ionization and harmonic generation in acetylene, using time-dependent density functional theory with various laser wavelengths.

Findings

Ionization is maximized when molecules align with laser polarization.

Excited states enhance high harmonic generation yield.

Photoelectron spectra are significantly affected by excited state populations.

Abstract

Ionization of acetylene by linearly-polarized, vacuum ultraviolet (VUV) laser pulses is modelled using time-dependent density functional theory. Several laser wavelengths are considered including one that produces direct ionization to the first excited cationic state while another excites the molecules to a Rydberg series incorporating an autoionizing state. We show that for the wavelengths and intensities considered, ionization is greatest whenever the molecule is aligned along the laser polarization direction. By considering high harmonic generation we show that populating excited states can lead to a large enhancement in the harmonic yield. Lastly, angularly-resolved photoelectron spectra are calculated which show how the energy profile of the emitted electrons significantly changes in the presence of these excited states.

Tables

Table 1. Table 1: Static properties of acetylene. The equilibrium C–C and C–H bond lengths and ionization potential calculated using LDA-PW92-ADSIC are compared with experimental values. The vertical ionization potentials from the X 1 ​ Σ g + superscript 𝑋 1 superscript subscript Σ 𝑔 X^{1}\Sigma_{g}^{+} ground state to the X 2 ​ Π u superscript 𝑋 2 subscript Π 𝑢 X^{2}\Pi_{u} , A 2 ​ Σ g + superscript 𝐴 2 superscript subscript Σ 𝑔 A^{2}\Sigma_{g}^{+} and B 2 ​ Σ u + superscript 𝐵 2 subscript superscript Σ 𝑢 B^{2}\Sigma^{+}_{u} cationic states are estimated, respectively, from the HOMO, HOMO-1 and HOMO-2 orbital energies. Experimental bond lengths are taken from Yasuike & Yabushita  Yasuike and Yabushita ( 2000 ) while the ionization potentials are taken from Wells & Lucchese  Wells and Lucchese ( 1999 ) .

	Equilibrium bond lengths (Bohr)		Vertical Ionization Potentials (Hartree)
	C–C bond length	C–H bond length	$X^2{\Pi }_u$	$A^2{\Sigma }_g^{+}$	$B^2{\Sigma }_u^{+}$
Current	2.207	2.045	0.4150	0.5550	0.6118
Experimental	2.273	2.003	0.4191	0.6140	0.6912

Equations

$$n(\mbox{\boldmath$r$},t)=2\sum_{j=1}^{N}\left|\psi_{j}(\mbox{\boldmath$r$},t)\right|^{2}.$$

$$\displaystyle i\frac{\partial}{\partial t}\psi_{j}(\mbox{\boldmath$r$},t)$$

$$\displaystyle=H_{\mbox{\scriptsize ks}}\psi_{j}(\mbox{\boldmath$r$},t)\hskip 56.9055ptj=1,\dots,N.$$

$$V_{\mbox{\scriptsize ext}}(\mbox{\boldmath$r$},\mbox{\boldmath$R$},t)=V_{\mbox{\scriptsize ions}}(\mbox{\boldmath$r$},\mbox{\boldmath$R$},t)+U_{\mbox{\scriptsize elec}}(\mbox{\boldmath$r$},t).$$

$$U_{\mbox{\scriptsize elec}}(\mbox{\boldmath$r$},t)=U_{L}(\mbox{\boldmath$r$},t)=\mbox{\boldmath$r$}\cdot\mbox{\boldmath$E$}(t),$$

$$U_{\mbox{\scriptsize elec}}(\mbox{\boldmath$r$},t)=U_{V}(\mbox{\boldmath$r$},t)=-i\mbox{\boldmath$A$}(t)\cdot\mbox{\boldmath$\nabla$}.$$

$$\mbox{\boldmath$A$}(t)=A_{0}\sin^{2}\left(\frac{\pi t}{T}\right)\cos(\omega_{L}t)\hat{{\mbox{\boldmath$e$}}}.$$

$$\mbox{\boldmath$E$}(t)=E(t)\mbox{\boldmath$\hat{e}$},$$

$$E(t)=E_{0}f(t)\sin ({\omega }_{L}t)-\frac{E_0}{{\omega }_L}\frac{\partial f}{\partial t}\cos ({\omega }_{L}t),$$

$$\displaystyle M_{k}\ddot{\mbox{\boldmath$R$}}_{k}=$$

$$\displaystyle-\frac{\partial}{\partial\mbox{\boldmath$R$}_{k}}\Bigr{(}V_{nn}(\mbox{\boldmath$R$})+Z_{k}\mbox{\boldmath$R$}_{k}\cdot\mbox{\boldmath$E$}(t)\Bigr{)},$$

$$\displaystyle\psi_{j}(\mbox{\boldmath$r$},t)$$

$$\displaystyle=\psi^{B}_{j}(\mbox{\boldmath$r$},t)+\psi^{S}_{j}(\mbox{\boldmath$r$},t).$$

$$M(\mbox{\boldmath$r$})=M_{x}(x)M_{y}(y)M_{z}(z).$$

$$M_x(x)=\{\begin{array}{cc}1& \left|x\right|\le x_m\\ 1-\alpha {\left(\left|x\right|-x_m\right)}^5& \left|x\right|>x_m\end{array},$$

$$\alpha =\frac{1-M_f}{(x_f-x_m{)}^5}.$$

$$S_{k}(\omega)=\left|\int e^{i\omega t}\,\hat{\mbox{\boldmath$e$}}_{k}\cdot\ddot{\mbox{\boldmath$d$}}(t)\,\mathrm{d}t\right|^{2},$$

$$\ddot{\mbox{\boldmath$d$}}(t)=-\int n(\mbox{\boldmath$r$},t)\mbox{\boldmath$\nabla$}H_{\mbox{\scriptsize ks}}\mathrm{d}\mbox{\boldmath$r$}.$$

$$\ddot{\mbox{\boldmath$d$}}_{j}(t)=-\int n_{j}(\mbox{\boldmath$r$},t)\mbox{\boldmath$\nabla$}H_{\mbox{\scriptsize ks}}\mathrm{d}\mbox{\boldmath$r$}=-2\int\left|\psi_{j}(\mbox{\boldmath$r$},t)\right|^{2}\mbox{\boldmath$\nabla$}H_{\mbox{\scriptsize ks}}\mathrm{d}\mbox{\boldmath$r$}.$$

$$\displaystyle b_{j}(\mbox{\boldmath$k$})\ e^{iT\mbox{\boldmath$k$}^{2}/2}$$

$$\displaystyle\equiv\int_{|\mbox{\boldmath$r$}|>R_{S}}d\mbox{\boldmath$r$}\ \chi_{\mbox{\boldmath$k$}}^{*}(\mbox{\boldmath$r$})\psi_{j}(\mbox{\boldmath$r$},T).$$

$$\mathcal{H}_{V}(t)=\frac{1}{2}\ \Biggr{[}\mbox{\boldmath$k$}-\mbox{\boldmath$A$}(t)\Biggr{]}^{2}.$$

$$\chi_{\mbox{\boldmath$k$}}(\mbox{\boldmath$r$})=\frac{1}{\left(\sqrt{2\pi}\right)^{3}}\ e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$r$}}e^{-i\Phi(\mbox{\boldmath$k$},t)},$$

$$\Phi(\mbox{\boldmath$k$},t)=\frac{1}{2}\int_{0}^{t}dt^{\prime}\ \Biggr{[}\mbox{\boldmath$k$}-\mbox{\boldmath$A$}(t^{\prime})\Biggr{]}^{2},$$

$$\displaystyle\left<\chi_{\mbox{\boldmath$k$}}\left|\theta(R_{S}){\left|\psi_{j}(\mbox{\boldmath$r$},T)\right.}\right.\right>$$

$$\displaystyle=i\int_{0}^{T}\mathrm{d}t\Bigr{<}\chi_{\mbox{\boldmath$k$}}\Bigr{|}\mathcal{H}_{V}(t)\theta(R_{S})-\theta(R_{S})\mathcal{H}(t)\Bigr{|}\psi_{j}(\mbox{\boldmath$r$},t)\Bigr{>}$$

$$\displaystyle=i\int_{0}^{T}dt\Biggr{<}\chi_{\mbox{\boldmath$k$}}\Biggr{|}\left[-\frac{1}{2}\nabla^{2}-i\mbox{\boldmath$A$}(t)\cdot\mbox{\boldmath$\nabla$},\,\theta(R_{S})\right]\Biggr{|}\psi_{j}(\mbox{\boldmath$r$},t)\Biggr{>},$$

$$e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$r$}}=4\pi\sum_{lm}i^{l}j_{l}(kr)Y_{lm}(\hat{\mbox{\boldmath$k$}})Y_{lm}^{\star}(\hat{\mbox{\boldmath$r$}}),$$

$$\displaystyle\Biggr{<}\chi_{\mbox{\boldmath$k$}}\Biggr{|}\left[-\frac{1}{2}\nabla^{2}-i\mbox{\boldmath$A$}(t)\cdot\mbox{\boldmath$\nabla$},\,\theta(R_{S})\right]\Biggr{|}\psi_{j}(\mbox{\boldmath$r$},t)\Biggr{>}=\sqrt{\frac{2}{\pi}}e^{i\Phi(\mbox{\boldmath$k$},t)}R_{S}^{2}\sum_{lm}(-i)^{l}Y_{lm}(\hat{\mbox{\boldmath$k$}})$$

$$\displaystyle\left[\begin{array}[]{l}\displaystyle-\frac{1}{2}\left(-\left.\frac{\partial j_{l}(kr)}{\partial r}\right|_{R_{S}}\psi_{j}(R_{S},t)+j_{l}(kR_{S})\left.\frac{\partial\psi_{j}(R_{S},t)}{\partial r}\right|_{R_{S}}\right)\\[22.76228pt] \displaystyle-iA_{x}j_{l}(kR_{S})\left<Y_{lm}(\hat{\mbox{\boldmath$r$}})\left|\sin\theta\cos\phi\left|\psi_{j}(R_{S},t)\right.\right.\right>\\[11.38092pt] \displaystyle-iA_{y}j_{l}(kR_{S})\left<Y_{lm}(\hat{\mbox{\boldmath$r$}})\left|\sin\theta\sin\phi\left|\psi_{j}(R_{S},t)\right.\right.\right>\\[11.38092pt] \displaystyle-iA_{z}j_{l}(kR_{S})\left<Y_{lm}(\hat{\mbox{\boldmath$r$}})\left|\cos\theta\left|\psi_{j}(R_{S},t)\right.\right.\right>\end{array}\right].$$

Full text

I Introduction

Attoscience is a field that has seen spectacular progress over the last several years with the continuing development of attosecond laser sources with which to probe dynamical processes in molecules, nanoscale devices and solids Krausz and Ivanov (2009). Since electrons move on the attosecond timescale, one of the main goals of attosecond chemistry is steering reactions through correlated electronic processes Lépine et al. (2014). Realizing this goal is key to the development of future ultrafast technologies: examples include the design of electronic devices Joachim et al. (2000), probes and sensors Willard and van Orden (2003), biological repair and signalling processes and development of optically-driven ultrafast electronics Krausz and Stockman (2014). Studying attosecond processes in molecules involves initiating a reaction by driving electrons far from equilibrium using a pump laser: one of the most widely used techniques is to ionize the molecule. The subsequent evolution of the electronic dynamics is followed by probing the system at later times using methods such as high-harmonic spectroscopy and photoelectron spectroscopy Lépine et al. (2014); Calegari et al. (2014); Trabattoni et al. (2015); Marangos (2016). Unravelling these dynamical processes presents a huge challenge due to the number of potential pathways that exist as well as the inherent difficulties of relating the initial and intermediate states to the photo-products produced. For this reason studying small molecules has attracted a lot of effort since they represent a compromise between having sufficient complexity and limiting the number of photo-products when compared to much larger molecules.

Acetylene (C2H2) is a small, linear polyatomic molecule which is isolelectronic to N2. It has the ground configuration $(1\sigma_{g})^{2}(1\sigma_{u})^{2}(2\sigma_{g})^{2}(2\sigma_{u})^{2}(3\sigma_{g})^{2}(1\pi_{u})^{4}$. In addition, the next lowest unoccupied orbitals are $(1\pi_{g})^{0}(3\sigma_{u})^{0}(4\sigma_{g})^{0}$. Ionization of acetylene has been studied both theoretically and experimentally by many groups. This is in part due to its rich electronic structure and in part due to its importance in many areas of chemistry. Initial studies of acetylene considered ionization by electron Collin and Delwiche (1967); King and Price (2007) and photon Collin and Delwiche (1967); Langhoff et al. (1981); Machado et al. (1982); Holland et al. (1983); Yasuike and Yabushita (2000); Fronzoni et al. (2003) impact. More recently, several studies have considered the response to strong fields. This has included single Argenti et al. (2012); Ji et al. (2015) and double Gaire et al. (2014); Gong et al. (2014) ionization, dissociative ionization of highly excited states Zammith et al. (2003), high harmonic generation Vozzi et al. (2010); Torres et al. (2010); Vozzi et al. (2012); Negro et al. (2014) and more recently initiating isomerization of acetylene to vinylidene Ibrahim et al. (2014).

The theoretical and computational study of molecules interacting with intense, ultrashort pulses is highly demanding. Several techniques have been developed to model laser-molecule interactions. These include simple yet powerful techniques like the strong field approximation (SFA) Lewenstein et al. (1994); Spanner and Patchkovskii (2009); Smirnova et al. (2009) and quantitative rescattering (QRS) theory Le et al. (2009); Lin et al. (2010) for modelling high-harmonic generation (HHG) in complex molecules. However, these approaches suffer from the drawback that only a small number of channels are generally included in a calculation. This can become an issue, depending on the molecule under investigation, if other channels become important during the interaction with the laser. In order to overcome these shortcomings, a range of techniques, based on ab initio methods, have been developed. Pioneering ab initio studies using density functional theory (DFT) Fronzoni et al. (2003); Remacle and Levine (2006); Argenti et al. (2012); Calegari et al. (2014) and configuration interaction (CI) methods Lünnermann et al. (2008) have been used to study electron migration in molecules following sudden removal of an electron and have observed charge oscillations along the full length of the molecule with a period of several femtoseconds. The drawback of these approaches is that the initial pumping of the molecule by a laser pulse is not described. In light of this several approaches that describe the explicit interaction with the laser field have been developed Thachuk et al. (1996); Mitrić et al. (2009); Richter et al. (2011).

One ab initio method that is widely used to treat electronic dynamics is time-dependent density functional theory (TDDFT) Runge and Gross (1984). The TDDFT method has been extensively applied to the study of laser-molecule and laser-cluster interactions Kunert and Schmidt (2003); Calvayrac et al. (2000); Marques et al. (2003). Since we generally want to include fragmentation processes in molecules we must couple the treatment of the electronic subsystem to another approach for treating the ionic dynamics. One popular approach is the Ehrenfest method whereby the quantum treatment of the electrons is coupled to a classical treatment of the ions, resulting in a method known as non-adiabatic quantum molecular dynamics (NAQMD) Kunert and Schmidt (2003).

In this paper we study ionization processes in acetylene using this mixed quantum-classical method. The paper is laid out as follows. In section II we describe our method and show how it is applied to the current study. In particular we describe how we calculate observables such as photoelectron spectra (PES). In section III we apply our method to study ionization of acetylene by VUV laser pulses. After obtaining the equilibrium ground state we firstly consider how the wavelength and orientation of the laser pulse alters ionization. Then we show that resonant excitation of particular excited states can dramatically enhance HHG. Lastly, we calculate angularly-resolved PES for acetylene at two distinct wavelengths and show how excited states significantly alter the emitted electrons. Finally, some conclusions and directions for future work are given in section IV.

II Theoretical Approach

In our calculations we implement the NAQMD method in a code called EDAMAME (Ehrenfest DynAMics on Adaptive MEshes) which is described in more detail in previous papers Dundas (2012); Wardlow and Dundas (2016). We now briefly describe how it is applied to the calculations considered here. In sections II.1 and II.2 we set out the details of the method for treating the electrons and ions respectively. Section II.3 gives information on the numerical implementation of the approach, highlighting many of the parameters used for the calculations presented. Section II.4 describes how we calculate observables such as the amount of ionization, the HHG spectra and angularly-resolved photoelectron spectra.

II.1 TDDFT description of the electronic dynamics

We consider a molecule consisting of $N_{e}$ electrons and $N_{n}$ classical ions where $M_{k}$, $Z_{k}$ and \mbox{\boldmathR}_{k} denote respectively the mass, charge and position of ion $k$. Additionally \mbox{\boldmathR}=\{\mbox{\boldmathR}_{1},\dots,\mbox{\boldmathR}_{N_{n}}\}. Neglecting spin effects, the time-dependent electron density for the electrons can be expressed in terms of $N=N_{e}/2$ time-dependent Kohn-Sham orbitals, \psi_{j}(\mbox{\boldmathr},t), each having an initial occupancy of two, as

[TABLE]

These Kohn-Sham orbitals satisfy the time-dependent Kohn-Sham equations (TDKS)

[TABLE]

In Eq. (2), V_{\mbox{\scriptsize H}}(\mbox{\boldmathr},t) is the Hartree potential, V_{\mbox{\scriptsize ext}}(\mbox{\boldmathr},\mbox{\boldmathR},t) is the external potential, and V_{\mbox{\scriptsize xc}}(\mbox{\boldmathr},t) is the exchange-correlation potential. Both the Hartree and exchange-correlation potentials are time-dependent due to their functional dependence on the time-dependent density. The external potential accounts for both electron-ion interactions and the interaction of the laser field with the electrons and can be written as

[TABLE]

For the laser pulses considered in this paper, the innermost electrons can be considered frozen and so we need only consider the response of the valence electrons. In that case, we describe the electron-ion interactions, V_{\mbox{\scriptsize ions}}(\mbox{\boldmathr},\mbox{\boldmathR},t), with Troullier-Martins pseudopotentials Troullier and Martins (1991) in the Kleinman-Bylander form Kleinman and Bylander (1982). The pseudopotentials were generated using the Atomic Pseudopotential Engine (APE) Oliveira and Nogueira (2008).

All calculations are carried out using the dipole approximation. Within EDAMAME, both length and velocity gauge descriptions of the electron-laser interaction, U_{\mbox{\scriptsize elec}}(\mbox{\boldmathr},t), can be used. For the length gauge we have

[TABLE]

while for the velocity gauge

[TABLE]

In these equations the vector potential, \mbox{\boldmathA}(t), is related to the electric field, \mbox{\boldmathE}(t), through the relation ${\bf E}(t)=-\frac{\partial}{\partial t}{\bf A}(t)$. Therefore, for the linearly polarised laser pulses used in this work we consider a $\sin^{2}$ pulse envelope and write

[TABLE]

Here, $A_{0}$ is the peak value of the vector potential, $\omega_{L}$ is the laser frequency, $T$ is the pulse duration and \hat{\mbox{\boldmathe}} is the unit vector in the polarization direction of the laser field. From this, the electric vector is

[TABLE]

where

[TABLE]

and where $E_{0}$ is the peak electric field strength. We find good agreement between the solutions using both gauges. However, using the length gauge is computationally less expensive than the velocity gauge and so in most calculations we will use the length gauge. We will use the velocity gauge for calculating the photoelectron spectra using the t-SURFF method (our implementation of t-SURFF in EDAMAME is described later).

Within TDDFT all electron correlation effects are included through the exchange-correlation potential, V_{\mbox{\scriptsize xc}}(\mbox{\boldmathr},t), which is derived from an exchange-correlation action functional. This functional is unknown and so must be approximated. In addition it contains information on the whole history of past densities. In order to avoid these memory effects the adiabatic approximation is generally used so that ground-state approximations to the exchange-correlation functional can be employed Burke et al. (2005). In our calculations we use the local density approximation (LDA) incorporating the Perdew-Wang parameterization of the correlation functional Perdew and Wang (1992). The LDA describes the electrons as a homogeneous electron gas. Although this functional is widely used, it suffers from self-interaction errors which means that it does not have the correct long-range behaviour. One major consequence is that electrons are too loosely bound and excited states are not accurately described Casida et al. (1998). Therefore, we supplement this functional with the average density self-interaction correction (ADSIC) Legrand et al. (2002) which reinstates the correct long-range behaviour in an orbital-independent fashion. This implementation will be referred to as LDA-PW92-ADSIC. The description of autoionizing resonances arising from double excitations is problematic in TDDFT. However, in the present work resonance phenomena primarily arise from single excitations which are described well using TDDFT Krueger and Maitra (2009). In particular, previous TDDFT studies of photoionization of acetylene have shown a good level of accuracy Fronzoni et al. (2003).

II.2 Classical description of the ionic dynamics

The ionic dynamics are treated classically using Newton’s equations of motion. For ion $k$ we have

[TABLE]

where V_{nn}(\mbox{\boldmathR}) is the Coulomb repulsion between the ions and Z_{k}\mbox{\boldmathR}_{k}\cdot\mbox{\boldmathE}(t) denotes the interaction between ion $k$ and the laser field. A major benefit of LDA-PW92-ADSIC is that the exchange-correlation potential is obtained as the functional derivative of an exchange-correlation functional, meaning that the forces acting on ions can be derived from the Hellmann-Feynman theorem.

II.3 Numerical implementation of the method

We solve the Kohn-Sham equations numerically using finite difference techniques on a three dimensional Cartesian grid Dundas (2012). All derivative operators are approximated using 9-point finite difference rules and the resulting grid is parallelized in 3D. Converged results were obtained using grid spacings of $0.4\,a_{0}$ for all coordinates. When using linearly-polarized laser pulses we have to treat ionizing and recolliding wavepackets that mainly travel along the laser polarization direction. We align the pulse along $z$ and use a larger grid extent in this direction.

Both the TDKS equations, Eq. (2), and ionic equations of motion, Eq. (9), need to be propagated in time. For the TDKS equations, this is achieved using an 18th-order unitary Arnoldi propagator Dundas (2004); Arnoldi (1951); Smyth et al. (1998). For the ionic equations of motion, we propagate using the velocity-Verlet method. Converged results are obtained for a timestep of $0.2$ a.u. for both sets of equations.

During the interaction of the molecule with the laser pulse, ionizing wavepackets can travel to the edge of the spatial grid, reflect and travel back towards the molecular centre. These unphysical reflections can be avoided using a wavefunction splitting technique that removes those wavepackets reaching the edges of the grid Smyth et al. (1998). The splitting is accomplished using a mask function, M(\mbox{\boldmathr}), that splits the Kohn-Sham orbital, \psi_{j}(\mbox{\boldmathr},t), into two parts

[TABLE]

In this equation, the first term is located near the molecule and is associated with non-ionized wavepackets. The second term is located far from the molecule and is associated with ionizing wavepackets: this part is discarded in our calculations. The point at which we apply this splitting must be be chosen carefully to ensure only ionizing wavepackets are removed. We write the mask function in the form

[TABLE]

Such a decomposition is most appropriate when a grid having different extents in each direction is used. If we consider the $x$ component we can write

[TABLE]

where

[TABLE]

Here $x_{m}$ is the point on the grid where the mask starts, $x_{f}$ is the maximum extent of the grid in $x$ and $M_{f}$ is the value that we want the mask function to take at the edges of the grid. Similar descriptions are used for $M_{y}(y)$ and $M_{z}(z)$.

II.4 Calculation of observables

Within TDDFT all quantities, including observables, are functionals of the electronic density. As in the case of the exchange-correlation functional, the exact form of these functionals are unknown in many cases and approximations must be made. In this paper, we present results for ionization, high-harmonic generation and photoelectron spectra. We describe below how these observables are calculated.

II.4.1 Ionization

The exact form of this functional is unknown and so most measures of ionization generally are obtained using geometric properties of the time-dependent Kohn-Sham orbitals Ullrich (2000). In this approach, bound- and continuum-states are separated into different regions of space through the introduction of an analysing box. In principal the continuum states occupy the regions of space where the wavefunction splitting is applied and hence the reduction in the orbital occupancies provides a measure of ionization.

II.4.2 High harmonic generation

For HHG, we calculate the spectral density, $S_{k}(\omega)$, along the direction \hat{\mbox{\boldmathe}}_{k} from the Fourier transform of the dipole acceleration Burnett et al. (1992)

[TABLE]

where \ddot{\mbox{\boldmathd}}(t) is the dipole acceleration given by

[TABLE]

Additional information can be gained by calculating the response for each Kohn-Sham orbital Chu and Groenenboom (2016). In that case we calculate the dipole acceleration for each state as

[TABLE]

While this neglects interferences between different orbitals in the overall harmonic signal, it does give an indication of the the contribution of that state.

II.4.3 Photoelectron Spectra

Photoelectron spectra are calculated using the time-dependent surface flux (t-SURFF) method. t-SURFF has been developed by Tao and Scrinzi for one and two electron systems Tao and Scrinzi (2012); Scrinzi (2012), and it has been further applied to the strong field ionization of H${}_{2}^{+}$ Yue and Madsen (2014, 2013), and more complex systems Karamatskou et al. (2014); Wopperer et al. (2017). In this method we think of configuration space as being divided into two regions, divided by a spherical surface at radius $R_{S}$. This is similar to the situation described by Eq. (10), with the splitting occurring at $R_{S}$. For a sufficiently long time, $T$, after the interaction with the laser pulse, we assume that the inner region only contains the bound state wavepackets whereas the outer region only contains the ionized wavepackets. The scattering amplitudes, b_{j}(\mbox{\boldmathk}), for a given Kohn Sham orbital, \psi_{j}(\mbox{\boldmathr},T), can then be obtained by projecting onto scattering solutions, \chi_{\mbox{\boldmathk}}. These amplitudes are expressed as

[TABLE]

The key to the method is to convert the volume integral on the right-hand-side of Eq. (17) into a time-dependent surface integral. In order to do that, we must know the evolution of the wavefunction after it has passed through the surface. The Hamiltonian which rules the dynamics in the outer region is the Volkov Hamiltonian, which has the form

[TABLE]

The solutions to the corresponding TDSE are the Volkov solutions

[TABLE]

where

[TABLE]

is the Volkov phase and \mbox{\boldmathA}(t) is the vector potential of the electric field.

In order to transform the volume integral in Eq. (17) to the time integral of the surface flux we write

[TABLE]

where we have introduced the Heaviside function, $\theta(x)$. We decompose the plane wave in the Volkov solutions in terms of spherical harmonics in $k$- and $r$-space as

[TABLE]

where Y_{lm}(\hat{\mbox{\boldmathx}}) is a spherical harmonic in the angular variables of \hat{\mbox{\boldmathx}}-space for orbital angular momentum quantum number $l$ and magnetic quantum number $m$ and $j_{l}(x)$ is a spherical Bessel function of the first kind. The integrand in Eq. (21) is then calculated according to the formula Zielinski (2016)

[TABLE]

To implement t-SURFF, we need to know the wavefunction and its first derivative at a series of points on a spherical shell at $R_{S}$. In light of the spherical nature of the problem and the fact that the initial implementation of t-SURFF was in spherical coordinates, this is the most appropriate coordinate system to use. However, EDAMAME is written is Cartesian coordinates and so we use tri-cubic interpolation Lekien and Marsden (2005) to interpolate the Kohn-Sham orbitals onto a spherical shell: this is carried out during a simulation, output to file and post-processed afterwards. To obtain the photoelectron spectra we integrate \left|b_{j}(\mbox{\boldmathk})\right|^{2} over all angular variables and sum up the contribution of all orbitals. We can also calculate angularly-resolved photoelectron spectra. In a spherical representation the angular variables of $k$ are $(\theta_{k},\phi_{k})$. We integrate over the $\phi_{k}$ angle to produce plots of the spectrum for $\theta_{k}$ only. We consider two cases: one where the integration in $\phi_{k}$ is carried out over the positive $x$-axis, the other where it is carried out over the negative $x$-axis. Each of these will give a semi-circular plot which we merge together to form a circular plot. This allows us to take account of asymmetries in each region. This implementation of t-SURFF is freely available n a library called POpSiCLE (PhOtoelectron SpeCtrum library for Laser-matter intEractions) POpSiCLE library .

III Results

In this section we apply our method to study ionization of acetylene with ultra-short linearly-polarized VUV laser pulses. In particular we will study resonance phenomena arising from transitions to excited states. The section is laid out as follows. In section III.1 we calculate the equilibrium ground state of acetylene using our TDDFT method. Starting from this initial state we investigate the role of excited states in the response of acetylene to intense laser pulses. Firstly, in section III.2 we show how ionization to the various cationic states of acetylene changes as a function of the laser wavelength and the orientation between the molecule and the laser pulse. Secondly, in section III.3 we show how using a VUV pulse tuned to a particular excited state can give rise to a large enhancement in the HHG efficiency. Finally, in section III.4 we investigate the role of this excited state by analysing the photoelectrons produced after interaction with VUV pulses. In particular, two wavelengths are considered, one of which is tuned to a particular excited state, the other which bypasses it. In the calculations presented here the ions are allowed to move. The effect of allowing the ions to move, when studying HHG in benzene, was considered in a previous paper Wardlow and Dundas (2016). For the results presented here, in which acetylene interacts with IR pulses, we see the hydrogen atoms oscillating with amplitudes of around 0.1 a0, whereas when the molecule interacts with VUV pulses the hydrogen atoms oscillate with amplitudes around an order of magnitude smaller.

III.1 Equilibrium ground-state properties

The equilibrium ground state configuration of acetylene was obtained using a geometry relaxation scheme. Using the LDA-PW92-ADSIC exchange-correlation potential, this gave the equilibrium bond lengths and ionization potentials as shown in Table 1. Using Koopman’s theorem, the vertical ionization potentials are estimated from the energies of the Kohn-Sham orbitals. We see that these quantities are in good agreement with the experimental values. While the C–C bond length is well reproduced, the C–H bond length is slightly overestimated. Likewise the ionization potential to the ground cationic state is well reproduced using LDA-PW92-ADSIC. However, those for the next two excited cationic states are lower than experimental values. In addition to ionization (bound-continuum transitions), an incident laser pulse can also excite an electron to a bound excited state (bound-bound transitions). This corresponds to a transition from an occupied state to an unoccupied state. The Kohn-Sham orbital energies of the five active occupied orbitals, together with those for the several lowest unoccupied orbitals, are plotted in Fig. 1. These energies are obtained through direct diagonalization of the field-free Kohn-Sham Hamiltonian.

III.2 Ionization to different cationic states

The likelihood of ionization to a particular cationic state depends greatly upon the symmetry of the states involved (and thus the orientation of the molecule with the laser pulse) and the energy required for ionization. For example, consider ionization of acetylene by high-intensity IR laser pulses. In that case ionization from the ground state of the neutral molecule ($X^{1}\Sigma_{g}^{+}$) to the ground state of the cation ($X^{2}\Pi_{u}$), i.e. removal of an electron from the HOMO orbital, is more likely when the laser pulse is aligned perpendicular to the molecular axis rather than when it is aligned along the molecular axis. Likewise, ionization to the $A^{2}\Sigma_{g}^{+}$ cationic state is associated with removal of an electron from the HOMO-1 orbital and is more likely when the laser pulse is aligned along the molecular axis. Fig. 2 presents the population loss from each Kohn-Sham state during interaction of acetylene with 8-cycle (2.67 fs) linearly polarized VUV laser pulses having a peak intensity of $I=$ $1.0\times 10^{12}\mbox{ W/cm}^{2}$: in total fourteen different wavelengths were considered. This population loss is calculated for two orientations of the molecule with the laser pulses: parallel and perpendicular. Population loss from each Kohn-Sham orbital is associated with vertical ionization to different cationic states. We see that more ionization occurs in parallel alignment. This is the opposite of the response we observe at IR wavelengths using higher laser intsities. In the parallel orientation, we see large increases in the amount of ionization from each Kohn-Sham orbital at the photon energies associated with the orbital energies. In the perpendicular case the response is different, especially for the HOMO orbitals. In that case we see enhancement in ionization for a range of wavelengths.

III.3 The role of excited states in HHG

Previous studies of photoionization in acetylene have shown that the $3\sigma_{g}\rightarrow 3\sigma_{u}$ transition is associated with an autoionizing state which leads to a feature in the photoionization spectrum around 0.4963 Ha in the photon energy Yasuike and Yabushita (2000); Wells and Lucchese (1999). Such a transition is associated with an excitation from the HOMO-1 to the LUMO+1. Using linearly-polarized laser pulses with the polarization aligned parallel to the molecular axis, we have recently shown that this excitation can lead to a significant enhancement of the plateau harmonics Mulholland and Dundas . In our calculations the HOMO-1 to LUMO+1 energy gap is 0.4445 Ha. We consider a pump-probe scheme in which acetylene interacts with an 8-cycle (2.67 fs) linearly polarized VUV laser pulse having a wavelength of $\lambda=$ 102 nm (photon energy $=0.4467$ Ha) and a peak intensity of $I=$ $1.0\times 10^{12}\mbox{ W/cm}^{2}$. For this pulse, the bandwidth is sufficient to excite the $2\sigma_{u}\rightarrow 4\sigma_{g}$ transition as well. Immediately after the pump pulse, the molecule interacts with a 5-cycle (24.2 fs) laser pulse having a wavelength of $\lambda=$ 1450 nm (photon energy $=$ 0.0314  Ha) and a peak intensity of $I=$ $1.0\times 10^{14}\mbox{ W/cm}^{2}$. For these calculations the grid extents were $|x|=|y|\leq 90.8$ a0 and $|z|\leq 146.8$ a0.

Fig. 3 presents the harmonic response for this pump-probe scheme. We see that a double plateau structure is present. The outer plateau arises due to ionization from and recombination back to the HOMO-1. The inner plateau is due to ionization from and recombination back to the LUMO+1. However, the cut-off energy for this inner plateau is lower than that predicted by the classical three-step model: we attribute this to suppression of ionization from the LUMO+1 due to transitions back to the HOMO-1 Mulholland and Dundas . These transitions from the LUMO+1 to the HOMO-1 give rise to a window of enhanced harmonics (harmonics 11–21). We can show that the inner plateau is due to the $3\sigma_{g}\rightarrow 3\sigma_{u}$ transition rather than the $2\sigma_{u}\rightarrow 4\sigma_{g}$ transition by calculating the contribution of each state to the harmonic response using Eq. (16). The contributions of the HOMO-1 and HOMO-2 to the harmonic response are also shown in Fig. 3. We see that the HOMO-1 contribution dominates by over an order of magnitude.

To show that resonant excitation is the cause for the HHG enhancement, we can consider the response of acetylene to VUV pump pulses having different photon energies. Fig. 4 presents results for four different wavelengths. In each case the laser intensity of the pump pulse is kept constant at $I=$ $1.0\times 10^{12}\mbox{ W/cm}^{2}$ while the number of pulse cycles is changed so that the total pulse length is roughly the same for all four pulses ($\sim$ 2.5 fs). From this figure we see clear evidence of resonant excitation around $\lambda=$ 102 nm. The results for $\lambda=$ 82 nm are interesting. For this laser wavelength we are able to directly ionize the molecule to the $A^{2}\Sigma_{g}^{+}$ excited cationic state. In that case, we would expect significant ionization from the $3\sigma_{g}$ orbital meaning that the population in the LUMO+1 (excited from the $3\sigma_{g}$) should decrease. This clearly shows up in the plot as a reduction in the intensity of the bound-bound harmonics (harmonics 11–21).

III.4 Photoelectron spectra using VUV laser pulses

The previous results clearly show how excited states affect HHG in acetylene. Since these excitations were initiated by the interaction with a VUV pump pulse we now consider how excited states affect ionization. We will calculate photoelectron spectra for acetylene ionized by linearly polarized VUV laser pulses. In particular, two wavelengths will be considered: $\lambda=$ 82 nm and $\lambda=$ 102 nm. The shorter wavelength provides enough energy to ionize directly from the HOMO-1 orbital with one photon while the longer drives the $3\sigma_{g}\rightarrow 3\sigma_{u}$ transition with one photon. For these calculations the grid extents were $|x|=|y|\leq 76.8$ a0 and $|z|\leq 90.8$ a0. The surface was placed at $R_{s}=$ 60 $a_{0}$ and the maximum angular momentum quantum number used in Eq. (27) was $l_{\max}=100$. The Kohn-Sham equations were propagated for 20 fs after the laser pulse finishes.

For the shorter wavelength, we consider the interaction of acetylene with an 8-cycle linearly polarized VUV laser pulse having a wavelength of $\lambda=$ 82 nm and a peak intensity of $I=$ $1.0\times 10^{12}\mbox{ W/cm}^{2}$. Fig. 5 presents photoelectron spectra for two orientations of the molecule with the field. In Fig. 5(a) the pulse is aligned parallel to the molecular axis while in Fig. 5(b) the pulse is aligned perpendicular to the molecular axis. For this laser intensity only a few photons will be absorbed. The results show two broad peaks associated with one- and two-photon absorption. The width of each peak is consistent with the the bandwidth of the pulse. As well as plotting the full spectra (black lines) we also plot the contribution of each Kohn-Sham orbital. In the parallel orientation the HOMO orbitals have negligible response. For the one-photon peak we see that the HOMO-1 dominates whereas for the two-photon peak both the HOMO-1 and HOMO-2 contribute: in the lower energy region the HOMO-2 dominates while for the higher energy region the HOMO-1 dominates. This arises due to the energy gap between these two orbitals (marked on the graph). In the perpendicular case, the response is quite different. Firstly, the overall probability density is an order of magnitude smaller than in the parallel orientation. This is to be expected from the results in Fig. 2. Secondly, one of the HOMO orbitals has a significant response while the HOMO-2 orbital has little contribution to the spectrum. Lastly, the HOMO-1 still contributes strongly to the spectrum. Due to the large energy gap between the HOMO and HOMO-1 one- and two-photon peaks are very broad. Indeed, especially for the one-photon peak, the demarcation between each orbital’s contribution can clearly be seen.

Fig. 6 presents the angularly-resolved spectrum for the parallel orientation while Fig. 7 presents the angularly-resolved spectrum for the perpendicular orientation. In each figure we plot the full spectrum and the contribution from each Kohn-Sham orbital. The axes refer to the energy parallel to the pulse polarization direction and that perpendicular to it. Consider the parallel alignment results in Fig. 6. We can clearly see the one- and two-photon contribution from each state. For the HOMO orbitals electrons in the one-photon peak are emitted at 45o to the laser polarization while they are emitted at 35o and 90o for the two-photon peak. Consider the HOMO-1 and HOMO-2 orbitals. For the one-photon peak the HOMO-1 electrons are emitted at 55o while for the HOMO-2 electrons are mainly along the laser polarization. For the two-photon peak electrons are emitted predominantly along the polarization direction. For the HOMO-1 we also see some ionization at 55o and 90o whereas for the HOMO-2 see some emission at 70o. Now consider the perpendicular alignment results in Fig. 7. In this case we see that the response of the HOMO orbitals dominate. The HOMO-b responds predominantly and we see that electrons are emitted along the pulse polarization direction for both one- and two-photon processes. For the one-photon peak, HOMO-1 electrons are preferentially emitted along the pulse direction while HOMO-2 electrons are emitted at 55o.

For the longer wavelength, consider the interaction of acetylene with an 8-cycle linearly polarized VUV laser pulse having a wavelength of $\lambda=$ 102 nm and a peak intensity of $I=$ $1.0\times 10^{12}\mbox{ W/cm}^{2}$. Fig. 8 presents photoelectron spectra for the parallel orientation. From our HHG results we know that excited states are important for this orientation. Fig. 8(a) presents the full spectrum and the contribution of each Kohn-Sham orbital over the full energy range considered. While the basic features of the spectra are similar to those observed at $\lambda=$ 82 nm, there are a number of important differences. Firstly, the spectral density is lower. This is to be expected as one-photon absorption results in excitation of the HOMO-1 and HOMO-2 orbitals. More importantly, if we consider the two-photon peak we see that is highly structured. Fig. 8(b) zooms in on this peak and presents the spectrum calculated on a much finer energy mesh. The structure is a signature of the excited states populated during the absorption of one photon before subsequent ionization by a second photon. We clearly see a wide range of Rydberg states populated as well as several resonance features. For the HOMO-1 contribution this is generally assigned to the $3\sigma_{g}\rightarrow np\sigma_{u}$ Rydberg series whereas for the HOMO-2 it is assigned to the $2\sigma_{u}\rightarrow ns\sigma_{g}$ Rydberg series. Fig. 9 presents the angularly-resolved spectrum. Again we see that the response is quite different to that observed in Figure 6. In particular, if we consider the contribution of the HOMO-1 orbital we see that the two-photon peak has little angular variation. For the HOMO-2 orbital, on the other hand, electrons are predominantly emitted along the laser polarization direction.

IV Conclusions

We have presented calculations of ionization of acetylene by VUV laser pulses using a TDDFT approach and showed the important role that excited state dynamics have in the response. In particular we showed that orientation effects have a dramatic effect on the ionization dynamics. In this case enhanced ionization was found to occur whenever the laser is aligned parallel to the molecular axis. This is opposite to the situation found whenever high-intensity infra-red (IR) laser pulses are considered. In addition, by resonantly tuning the laser wavelength to these excited states, ionization and harmonic generation are significantly modified. In a recent paper we showed that HHG from acetylene can be enhanced whenever the molecule is excited to an autoionizing state Mulholland and Dundas . Here, we show that this behaviour is due to resonant excitation by considering a range of VUV wavelengths. Lastly, we have calculated angularly-resolved PES and showed how the spectrum of emitted electrons changes as the wavelength is tuned to various excitations of the molecule.

The excited states probed in this work incorporate resonances and autoionizing states. Such features are common to a wide range of molecules and should therefore alter the dynamical response to intense laser pulses. The method used in this work is of sufficient generality to be applied to other molecules: this includes more complicated molecules having fewer symmetries. As well as being able to study ionization and harmonic generation processes, the implementation of the t-SURFF method provides rich information on the energy distribution of the ionized electrons. A particular avenue for future work is studying circular dichroism in chiral molecules.

V Acknowledgments

AW and PM acknowledges financial support through PhD studentships funded by the UK Engineering and Physical Sciences Research Council. DD and AdlC acknowledges financial support from the European Union Initial Training Network CORINF and the embedded CSE programme of the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). This work used the ARCHER UK National Supercomputing Service and has been supported by COST Action CM1204 (XLIC).