Angular dependence of strong field sequential double ionization for neon and acetylene simulated with time-dependent configuration interaction using CIS and CISD-IP

Double ionization by strong laser fields and the correlated nature of the two-electron ejections have attracted much attention in the past two decades. These processes can be roughly categorized as sequential (SDI) and nonsequential (NSDI) depending on whether electron recollision plays a role.^1–3 However, such a delineation becomes problematic when the time of electron ejection becomes increasingly shorter toward less than 1 fsec. In such a process, the inherent electron correlation in the absence of electron recollision becomes important and should be treated explicitly. Therefore, modeling beyond the single active electron (SAE) approximation is critical for understanding multielectron effect in strong field double ionization of atomic and molecular systems. The recent experimental results suggest that ionization delays as short as 500 as can be accessed using the two-electron angular streaking method, in which the angle-dependent ionization rates were measured for both electrons to extract the time-resolved dynamics.^4,5 However, it is still difficult to extract electron correlation in this process with conventional methods.

To model such dynamics, the angular dependence of strong field single ionization has been simulated with rt-TDDFT^6–8 and TDCI approaches.^4,5,9–17 The former uses the time-dependent Kohn–Sham equations to propagate the electron density in real time in the presence of the oscillating field of the laser pulse. The latter uses the time-dependent Schrodinger equation to propagate the wavefunction which is expanded in terms of the field-free electronic configurations of the molecule. Complex absorbing potentials (CAP)¹⁸ are used in both approaches to absorb the density or wavefunction once it is far enough from the molecule and thereby simulate ionization. For single ionization, the simplest TDCI approach requires configuration interaction with all single excitations (CIS) and needs to include enough diffuse configurations to model the dynamics of the electron as it is propagated toward the complex absorbing boundary. We have used TDCIS with a complex absorbing potential to study the angular dependence of single ionization for a variety of systems.^4,5,9–17 Typical simulations involve propagating wavefunctions with thousands of singly excited configurations for tens of thousands of time steps.

Though accurate methods are available to simulate double ionization processes for two-electron systems,^19–40 these methods cannot be readily extended to multi-electron systems. Direct simulation of double ionization with the TDCI approach requires propagating a wavefunction with at least single and double excited configurations (CISD), which is practical only for small systems. For more typical molecules, millions of singly and doubly excited configurations would need to be propagated for tens of thousands of time steps. Since this is usually not practical, alternative methods need to be explored.

Sequential double ionization can be modeled by two separate TDCI or rt-TDDFT simulations. One electron is removed from the HOMO, HOMO−1, etc., leading to one or more states of the cation and yielding the probability of generating each cation state as a function of the laser field direction. Additional simulations for each of the relevant cation states are then used to model the different channels for the second ionization. We have used this approach to model the angular dependence of double ionization in benzene and in methyl iodide.^4,5 If the dynamics between close-lying or degenerate states of the cations is important, then these states should be treated together in one simulation, rather than in separate simulations. While the TDCIS-CAP approach is appropriate for the first step of a sequential double ionization, a more flexible representation may be desirable for simulating the second ionization step.

The CISD-ionization potential (CISD-IP) method described by Krylov and co-workers⁴¹ provides a means of representing a collection of cation states using a common set of orbitals and treating these cation states on an equal footing. The CISD-IP wavefunction consists of a small subset of the single and double excited configurations of the neutral molecule. The cation states are represented as ionizing excitations from the neutral ground state. These are augmented by a full set of singly excited configurations for each cation state (corresponding to double excitations from the neutral ground state) to allow for orbital relaxation. These excited configurations of the cations can be used in time-dependent CISD-IP with a complex absorbing potential (TDCISDIP-CAP) calculations to simulate the dynamics and ionization of the cations, analogous to the way that CIS provides the configurations needed for single ionization of the neutral molecule.

Experimental investigations of sequential double ionization of neon find that the second electron is preferentially ejected in a direction perpendicular to that of the first electron.^42,43 The measured angle-dependent ionization yields have been previously compared to model calculations for noble gas atoms.^44–46 Experimental studies of double ionization of aligned acetylene have determined angle-dependent ionization rates and fragment branching ratios.^47,48 Time-dependent density functional theory (TDDFT) calculations of the angle dependence of the first and second ionizations compared well with the experimental results.⁴⁸ As proof of concept, the present work uses the TDCIS-CAP and TDCISDIP-CAP methods to simulate the angular dependence for the second step of sequential double ionization of neon and acetylene, i.e., for the ionizing neon cation and acetylene cation. Section II summarizes simulations with the TDCI approach, describes the complex absorbing potential and additional basis functions used for modeling ionization, and outlines the calculation and analysis of the ionization rates. To model the first ionization of neon by linearly polarized light, we chose to remove an electron from the $2pyβ$ orbital. To avoid any spurious dynamics due to non-stationary states, the TDCI simulation of the second ionization was started from the eigenstate of the cation that is dominated by the $2pyβ$ hole. Ionization of the cation was found to occur primarily from a p orbital perpendicular to the p orbital involved in the first ionization, in agreement with the experimental observations. Similarly, for acetylene, the first electron was removed from the $πyβ$ orbital and the TDCI simulation of the second ionization was started from the eigenstate of the cation that is dominated by the $πyβ$ hole. In agreement with the experiment, ionization of the acetylene cation occurred primarily from the π orbital perpendicular to the π orbital involved in the first ionization.

In the TDCI-CAP approach, the time-dependent wavefunction is expanded on the basis of ground and excited configurations,

and is numerically integrated with the time-dependent Schrodinger equation,

Here, $H^el$ is the field-free molecular Hamiltonian for a fixed nuclear configuration. Molecular interaction with the field is described within the semiclassical dipole approximation $μ→^⋅E→τ$⁠, where $μ→^$ is the dipole operator and $E→τ$ is the time-dependent applied electric field. Ionization is captured with the complex absorbing potential (CAP), $iV^abs$⁠, which absorbs the wavefunction as an electron is propagated away from the molecule. Using the modified-midpoint approximation, the wavefunction is incrementally propagated with a timestep of Δt,

and the exponential is approximated with Trotter–Suzuki factorization.⁴⁹ The final equation of motion can be expressed in matrix notation as

where C(t) is a column vector of CI coefficients and d is a diagonal matrix of dipole eigenvalues obtained from d = WμW^T, where W is the matrix of eigenvectors.

The absorbing potential is constructed from a set of overlapping potentials on each atom as described in previous studies.^9–17 The value of the molecular absorbing potential is set to the minimum of values of the atomic absorbing potentials. The atomic potential for each atom is spherical and is zero for radial distances less than R₀, where R₀ is set to be 3.5 times the atomic van der Waals radius. The potential beyond R₀ rises quadratically to 5 hartree at approximately R₀ + 14.2 bohrs. The potential turns over quadratically to a constant value of 10 hartree at approximately R + 28.3 bohrs. In this work, R₀ is taken to be 10.72 bohrs for neon, 12.735 bohrs for carbon, and 9.544 bohrs for hydrogen (i.e., 3.5 times the van der Waals radius).

For simulations with TDCI-CAP, additional diffuse basis functions are required in addition to the standard basis set in order for the system to interact with the complex absorbing potential. Appropriate absorbing basis sets were developed and tested in previous studies.^10,14 For all simulations in this work, this absorbing basis set consists of four s functions with exponents of 0.0256, 0.0128, 0.0064, and 0.0032; four p functions with exponents of 0.0256, 0.0128, 0.0064, and 0.0032; four d functions with exponents of 0.0512, 0.0256, 0.0128, and 0.0064; and two f functions with exponents of 0.0256 and 0.0128. These basis functions were placed on each atom in addition to the standard aug-cc-pVTZ basis set.⁵⁰ 

The TDCI-CAP simulations for the second ionization are carried out using either the CIS or the CISD-IP states in Eq. (1). The TDCIS-CAP method employs the molecular orbitals and energies of the Hartree–Fock (HF) ground state of the cation [Scheme 1(a)]. The wavefunction is expanded with the corresponding ground HF and all singly excited determinants (CIS),

The CIS Hamiltonian matrix elements are given by

where E₀ is the ground state energy of the cation, ϵ are orbital energies, {i,j} refer to occupied molecular orbitals, and {a, b} refer to virtual molecular orbitals. The double bar integrals are

The one-electron matrix elements for the absorbing potential and the dipole are given by

where the one-electron integrals are $rA^s=ϕrA^ϕs=∫dr1ϕr*r1A^ϕsr1$⁠.

The CISD-IP Hamiltonian, as described by Krylov and co-workers,⁴¹ is constructed from the molecular orbitals and energies computed from the Hartree–Fock (HF) ground state of the neutral molecule, and includes all unique 1-hole, 2-hole, 1-particle determinants [Scheme 1(b)],

where the 1-hole determinants represent ionized configurations and 2-hole and 1-particle determinants represent singly excited ionized configurations. With these states, the matrix elements of the molecular Hamiltonian can be expressed as

where E₀ is the HF energy of the neutral ground state, {x,y} refers to ionized (hole) beta molecular orbitals, and i < x and j < y for the beta orbitals. The one-electron matrix elements are given by

In this work, the instantaneous ionization rate is defined as the change in the norm squared of the wavefunction,

The total ionization rate for a normalized CIS wavefunction in terms of the one-electron integral of the absorbing potential is

where $V00abs=∑iocciV^absi$⁠. The total rate can be partitioned into contributions to ionization from each of the occupied orbitals,

This can be interpreted as the rate of formation of an ionized state with a hole in orbital i.

For a normalized CISD-IP wavefunction, the total ionization rate is

This can be partitioned into the rate of formation of dication species with holes in orbitals x and i,

For field-free simulations, there is a small residual rate arising from the interaction of the ground state with the absorbing potential. This residual rate is subtracted from the rate in the field to obtain the net mean ionization rate,

Molecular orbitals and energies, one and two-electron integrals, and absorbing potential integrals were calculated with a locally modified copy of the Gaussian program suite.⁵¹ Closed shell systems were calculated with spin-restricted methods and open shell systems were calculated with spin-unrestricted methods. A stand-alone Fortran 95 code used in our previous TDCIS-CAP studies^4,5,9–17 was extended to include the CISD-IP Hamiltonian.⁴¹ This code was used to integrate the TDCI equations and analyze the TDCI results. The wavefunction was propagated in the oscillating electric field of a laser pulse for 16 000 steps using a time step of t = 0.05 a.u. (1.2 as), for a total simulation time of 19.35 fs (800 a.u.). The test in prior studies showed that reducing the time step by half changed the norm at the end of the simulation by less than 0.01%.^14,16 A 800 nm laser pulse was modeled by a linearly polarized oscillating electric field with a trapezoidal pulse,

with a period of $τp=2π/ω$ and maximum field strength of E_max was applied for 7 cycles for all systems. In this work, the mean ionization rates used in Eq. (17) correspond to the average of the instantaneous rates computed every 50 timesteps (60 as) for the duration of the sixth cycle of the field. For neon, the mean ionization rates in the z direction are 0.069 75 fs⁻¹ for TDCIS-CAP at E_max = 0.40 a.u. and 0.070 78 fs⁻¹ for TDCISDIP-CAP at E_max = 0.50 a.u. The norms squared after the pulse are 0.4133 and 0.4078, respectively; the residual field-free rates are 0.001 699 fs⁻¹ and 0.001 215 fs⁻¹, respectively. For C₂H₂, the mean ionization rates in the x direction are 0.093 92 fs⁻¹ for TDCIS-CAP at E_max = 0.105 a.u. and 0.045 54 fs⁻¹ for TDCISDIP-CAP at E_max = 0.115 a.u. The norms squared after the pulse are 0.3279 and 0.5748, respectively; the residual field-free rates are 0.000 796 8 fs⁻¹ and 0.000 965 11 fs⁻¹, respectively.

The electron dynamics of the ionization of the neon cation by a strong laser field were simulated with the TDCIS-CAP and TDCISDIP-CAP approaches. For the TDCIS-CAP simulation, the matrix elements of the Hamiltonian and the field-free singly excited states were computed from the molecular orbitals, energies, and integrals from a spin-unrestricted Hartree–Fock calculation of the cation with the hole localized in the $2pyβ$ orbital. For the TDCISDIP-CAP simulation, the Hamiltonian and field-free excited states were computed from a spin-restricted Hartree–Fock calculation of the closed shell neutral atom. The lowest three states of the CISD-IP Hamiltonian are degenerate and are dominated by a hole in the 2p_x, 2p_y, and 2p_z orbitals. The initial wavefunction for the TDCISDIP-CAP was set to be $Ψ0≈0.971Ψ2pyβ$⁠. In total, 869 CIS states and 3365 CISD-IP states were included in the simulations.

In Table I, the first ionization potentials (IP) calculated by ΔCCSD(T), the difference in the coupled cluster singles and doubles with perturbative triples [CISD(T)] energies of neon neutral and the neon cation are in very good agreement with the experimental values.⁵² The IP calculated by ΔSCF is lower than the IP calculated by the Koopmans theorem due to orbital relaxation in the SCF calculation of the cation. Since the TDCIS simulation starts with a Hartree–Fock calculation of the cation, the first IP listed for CIS is the same as the ΔSCF value. The IP calculated by CISD-IP is lower than the ΔSCF value because it includes some electron correlation for the cation. For second ionization potentials, there are three states to consider for the neon dication. Removing two electrons from the 2p^β orbitals of the neutral atom yields a component of the ³P state. Components of the ¹D and ¹S states can be obtained by removing an α electron and a β electron from the set of 2p orbitals. Spin-unrestricted SCF and CCSD(T) calculations of the triplet dication yield estimates for the IP from the cation to the ³P dication, while singlet spin-restricted SCF and CCSD(T) calculations of the dication yield IP estimates for the ¹D state but not the ¹S state. Koopmans theorem applied to the cation provides estimates for the second IP. Similar to the first IP, ΔCCSD(T) calculations for the second IP are in very good agreement with the experimental values for the ³P and ¹D dications.^53,54

Estimates of the second IP for the CIS and CISD-IP calculations can be obtained with the help of the absorbing potential matrix elements. The large absorbing basis set yields series of singly excited states that interact more strongly with the absorbing potential as they approach the ionization limit. Figure 1 shows the energies of the CIS and CISD-IP states and their absorbing potential matrix elements, $Vabs$⁠, as a function of the state number. The ionization potential for a particular series is obtained from the state energy where the absorbing potential reaches a maximum. The three series in Fig. 1 correspond to ionization of the cation to the ³P, ¹D, and ¹S dications. For CIS, the second IPs are about 1.5 eV lower than the Koopmans theory values, indicating that the diffused states are strongly absorbed a little before they reach the ionization limit. The second IPs for CISD-IP are about 10 eV higher than the CIS values for a combination of reasons. The CIS wavefunction uses the relaxed orbitals from the SCF calculation of the cation, whereas the CISD-IP wavefunction uses the unrelaxed orbitals from the neutral atom. In addition, the CISD-IP includes some electron correlation for the cation but not the dication. Since the second IPs are too high with CISD-IP, the calculated ionization rates will be too low. Therefore, we restrict the discussion to relative rates and the angular dependence of the ionization.

Figures 2(a) and 3(a) show the angular dependence of the second ionization of the $2pyβ$ neon cation in the yz plane (yx plane is equivalent by symmetry). The field strengths vary from 0.1 a.u. to 0.4 a.u. for TDCIS-CAP and 0.1 a.u. to 0.5 a.u. for TDCISDIP-CAP. For a given field strength, TDCISD-IP ionization rates are much smaller than those of TDCIS-CAP because of the difference in the second ionization potential calculated by these methods. However, the angular dependence of the ionization rate is similar, with a maximum when the field is polarized in the z direction and a minimum when the field is polarized in the y direction. This compares well with the angular dependence of the ionization rate extrapolated from the experimental data.⁴² 

The total ionization rate can be decomposed into contributions from the α and β orbitals using Eqs. (14) and (16) as shown in Figs. 2(b) and 2(c) for TDCIS-CAP and Figs. 3(b) and 3(c) for TDCISDIP-CAP. The initial state of the neon cation has an empty $2pyβ$ orbital. For both TDCIS-CAP and TDCISDIP-CAP, the largest contribution in the yz plane is ionization from the $2pzβ$ orbital yielding a component of the ³P dication state [Figs. 2(c) and 3(c)]. The next largest contribution in the yz plane is from ionization of the $2pzα$ orbital yielding a component of the ¹D dication state [Figs. 2(b) and 3(b)]. Corresponding contributions to the ³P and ¹D cations come from the $2pxβ$ and $2pxα$ orbitals in the xy plane. Ionization from the $2pyα$ orbital yields the ¹S dication. This component is larger in the TDCISDIP-CAP simulation and partially accounts for the smaller difference in second ionization rates between the y direction and z direction when compared to TDCIS-CAP.

Because of the larger number of configurations, the TDCISDIP-CAP simulation has some features that are lacking in the TDCIS-CAP simulation. Figure 4 shows the population of the 2s^β, $2pxβ$⁠, $2pyβ$⁠, and $2pzβ$ orbitals as a function of the field direction for five different times. Initially, the population of the $2pyβ$ orbital is near zero and the populations of the 2s^β, $2pxβ$⁠, and $2pzβ$ orbitals are near one. The population of the 2s^β and $2pxβ$ orbitals remains constant and independent of the angle of the field in the yz plane. For the field in the y or z directions, the population of the $2pyβ$ orbital stays small and the population of the $2pzβ$ orbital stays near one. However, for intermediate angles, the population of the $2pzβ$ orbital decreases with time and the population of the $2pyβ$ orbital increases. This corresponds to a rotation between degenerate states of the neon cation (i.e., change in the direction of the $2pβ$ hole) in response to the angle of the laser field. The polarizability is less in the direction of the $2pβ$ hole and greater in the direction perpendicular to it. The rotational force on the electron density or on the direction of the $2pβ$ hole resulting from the laser field is maximal at 45° to the axis of the $2pβ$ hole, similar to the rotational force on a homonuclear diatomic molecule in an electric field. CISDIP treats the interaction of the degenerate cation states with the field on an equal footing, but CIS treats only one of the degenerate states of the cation with the same accuracy as CISDIP.

Three dimensional views of the ionization yield are shown in Fig. 5. For a TDCIS-CAP simulation of the first ionization of an electron from the $2pyβ$ orbital of neutral neon, the angular dependence resembles the $2pyβ$ orbital as expected. In simulations of the second ionization of neon, the electron is ionized predominantly from orbitals perpendicular to the one involved in the first ionization, in this case the $2pxβ$ and $2pzβ$ orbitals. The TDCISDIP-CAP shows a bit more ionization in the y direction resulting from greater contributions from the α orbitals, as can be seen by comparing Figs. 2(b) and 3(b).

The ionization potentials for acetylene are listed in Table II. The lowest energy states of the cation and dication involve ionization of electrons from the π orbitals. The first IPs calculated by ΔCCSD(T) and Koopmans theorem are in very good agreement with the experimental value,⁵² while the ΔSCF and CISD-IP estimates are too low. The π² dication of acetylene is isoelectronic with O₂ and has ³Σ, ¹Δ, and ¹Σ as the three lowest states. Removing one electron from the $πxβ$ orbital and one electron from the $πyβ$ orbital of the neutral acetylene yields a component of the ³Σ state. Components of the ¹Δ and ¹Σ states can be obtained by removing one electron from a π^α orbital and one electron from a π^β orbital. The ΔCCSD(T) and Koopmans theory values for the second IP are in good agreement with the experimental value.^55–57 Similar to neon, the estimates of the second IP for CIS and CISD-IP calculations can be obtained by comparing the energies and the absorbing potential matrix elements, as shown in Fig. 6. The second IPs are about 1.5–2 eV too low for CIS and about 2 eV too high for CISD-IP.

The initial wavefunction for TDCIS-CAP started with the hole localized in the $πyβ$ orbital and TDCISDIP-CAP started from the state $Ψ0≈0.969Ψπyβ$⁠. The molecular axis of acetylene is placed along the z axis and the angular dependence of second ionization obtained from TDCIS-CAP and TDCISDIP-CAP simulations are shown in Figs. 7(a) and 8(a), respectively, for the xy plane. The TDCIS-CAP simulation used 2426 configurations and 7 field strengths varying uniformly from 0.085 a.u. to 0.115 a.u.; TDCSDIP-CAP simulations used 9420 configurations and the field strength from 0.095 a.u. to 0.125 a.u. Similar to the neon cation, the computed TDCIS-CAP ionization rates for the acetylene cation are larger than those of TDCISDIP-CAP due to the differences in the second IPs. The initial state of the acetylene cation was chosen to have an empty $πyβ$ orbital. For both the TDCIS-CAP and TDCISDIP-CAP simulations for ionization of the acetylene cation, the largest contribution in the xy plane is the removal of an electron from the $πxβ$ orbital yielding a component of the ³Σ state [Figs. 7(c) and 8(c)]. The next largest contribution in the xy plane is from ionization of the $πxα$ orbital yielding a component of the ¹Δ dication state [Figs. 7(b) and 8(b)]. Removing an electron from the $πyα$ yields a component of the ¹Σ dication.

Similar to the neon cation, the TDCISDIP-CAP simulation of ionization of the acetylene cation shows some evidence of rotation of the π^β hole. Figure 9 shows the population of the $πxβ$ and $πyβ$ orbitals as a function of the field direction for five different times. Initially the population of the $πxβ$ and $πyβ$ orbitals are approximately 1 and 0, respectively. For the field in the x or y direction, the populations of these orbital change very little. For intermediate angles, the population of the $πxβ$ orbital decreases with time and the population of the $πyβ$ orbital increases, corresponding to a rotation of the π ^β hole. However, the response is much smaller than that seen in the ionization of the neon cation (Fig. 4).

The angular dependence of the ionization yields for acetylene is shown in Fig. 10. The shape of the ionization yield for removing an electron from the $πyβ$ orbital resembles the shape of that orbital as expected. In the next ionization of this cation, the electron is removed predominantly from the $πxβ$ orbital which is perpendicular to the orbital involved in the first ionization. The TDCISDIP-CAP ionization yield is a little larger in the y direction, while the TDCIS-CAP ionization yield is a little larger in directions intermediate between x and z. Both show comparatively little ionization in the z direction, i.e., along the molecular axis, in keeping with the fact that the σ orbitals are significantly lower in energy than the π orbitals.

This study has demonstrated the use of TDCIS-CAP and TDCISDIP-CAP to simulate the second ionization of neon and acetylene. The estimated CIS second ionization potential is generally lower than CISD-IP because the CIS calculation starts with the relaxed orbitals of the cation, whereas the CISD-IP calculation starts with the orbitals of the neutral molecule. Since the second ionization potentials are larger for CISD-IP, higher intensities are needed to obtain second ionization rates comparable to the TDCIS-CAP simulations. Nevertheless, the two approaches yield similar angular dependence for the second ionization and show the second electron ionizing from the orbital perpendicular to the orbital involved in the first ionization. The TDCISDIP-CAP simulations of the angular dependence of the second ionization of neon and acetylene have some features arising from population (hole) dynamics not present in the TDCIS-CAP simulations. Because the TDCISDIP-CAP simulations include all single excitations from each of the degenerate monocation states, the field can induce a rotation between degenerate states of the cation. By contrast, the TDCIS-CAP simulations include single excitations from only one of the degenerate cation states and cannot model this effect. This new method puts us one step closer to modeling fully correlated strong field double ionization processes. In future work, we plan to link the simulations for the first and second ionizations so that the CI vector for the ionized state from a TDCIS simulation of the first ionization will be used as input for the TDCISDIP simulation for the second ionization.

This work was supported by a grant from the National Science Foundation (Grant No. CHE1856437). We thank the Wayne State University Computing Grid for access to high performance computing.