Simulation of photo-electron spectrum and electron scattering by dual time propagation

Investigation of electron and nuclear dynamics in intense laser pulses^1–22 has been an important research direction in the last two decades. Computer simulations describing these processes are very complicated because the laser steers the electrons far away from the atom and the electrons back-scatter on the ion later. The description of this process requires extremely large simulation cells in three dimensions. Examples for problems where large computational cells are required include (i) electron scattering, e.g., imaging molecules with laser-induced electron diffraction,^23–26 probing electron dynamics with attosecond photoelectron wave packets,^27,28 photon-induced near-field electron microscopy,²⁹ and (ii) photoelecton spectroscopy, e.g., photoexcitation circular dichroism or^19,30,31 photoelectron imaging.^32–34 

In a previous work, we have used the Volkov states³⁵ as the basis to solve the time-dependent Schrödinger-equation (TDSE)^36,37 for ionization problems. The approach is very simple and efficient, and one only has to propagate the potential energy matrix elements and the propagation is a simple multiplication with the time-dependent Volkov phase. The extension of the approach to large systems is difficult, however, because one has to calculate the potential matrix elements at every time step. In the present work, we bypass this problem by dividing the computational region into two parts.

In this work, we introduce a dual propagation approach where a central region will be described by the full time-dependent Hamiltonian and the outer region, where the laser field dominates, will be described by the Hamiltonian comprising the kinetic energy and the laser potential. The inner and outer regions are coupled by a source potential, and the electrons can enter or leave the central region without reflection. In the outer region, one can use Volkov states³⁵ to propagate the wave function. This propagation is numerically inexpensive, and the only limit on the size of the outer region is the available memory to store the coefficients of the wave function.

We will implement the dual propagation approach to the framework of the time-dependent density functional theory (TDDFT)³⁸ Hamiltonian. The TDDFT Kohn–Sham equations are often solved by time propagation of the ground state orbitals. Many approaches have been developed and tested for efficient time propagation.^39–48 

Considering the importance of the problem, it is not surprising that much research has been devoted to developing efficient computational methods for calculations in large simulation cells. The two most popular approaches are the complex absorbing potentials (CAPs) and the mask method.^49,50 In the CAP approach, an imaginary potential is added near the boundary acting as electronic density sinks.^51,52 The wave function that passes into regions with a CAP is lost, and the computational domain has to be made large to accommodate the CAP in realistic calculations. There are various forms of CAPs that are developed^53–55 to minimize the reflections.

In the mask method or wave function splitting technique, the system is divided into two regions: (A) a region treated quantum mechanically, usually via density functional theory (DFT), and (B) a region far from any nuclei, in which the Coulombic interactions with the nuclei and other electrons are negligible.^10,56–58 This splitting of domains is advantageous because in the absence of electron–nuclear and electron–electron interactions, the dynamics of a particle in an external time-dependent field may be described analytically using Volkov states³⁵ or free propagators.^49,50

The mask function method is similar to the dual propagation in the sense that both approaches divide the space into inner and outer regions and uses the Volkov states in the outer part. There are two problems with the mask method, which restricts its accuracy and applicability. The first is that the mask function is not perfect and spurious reflections may appear.⁵⁸ Second, to take care of the electrons reentering into the inner region, the momentum space wave functions have to be transformed back to real space. This step can only be done by discrete Fourier transformations, which introduce unwanted periodic boundary conditions.⁵⁸ This is highly undesirable because the electron density leaving the cell re-enters on the opposite side. These effects are especially problematic for stronger laser pulses. The problems of the mask method are discussed in details in Appendix B of Ref. 58 and in Ref. 59.

In the dual approach, on the other hand, the source potential connects the inner and outer wave functions, and there are no artificial reflections or periodicity at the boundary. The initial inner wave function is set to zero in the entire inner region at each time when it approaches the boundary region, and the source potential takes into account the electrons leaving and re-entering into the box.

We will test the approach calculating photoelectron spectrum (PES) after ionization and studying electron scattering on molecules. The photoelectron spectrum has been calculated by many different methods, including the surface flux method^18,60–62 where the time-dependent outgoing flux is monitored and the mask method.^63,64 These approaches have been tested, and their relative merits have been discussed in Refs. 65 and 66.

As an illustrative application, we will calculate the three-dimensional photoelectron angular distribution (PAD) of CO and C₃H₆O (metyloxirene) molecules in circularly polarized laser field. These calculations require very large computational cells due to the complex ionization dynamics. In the calculations, we use a selected laser-molecule orientation. To compare the simulations with experiments, one has to repeat the calculation with many different orientations as we did in Ref. 67 for hydrocarbon molecules in linearly polarized light.

Many physical interaction can lead to photoelecton emission ranging from single-photon ionization by extreme ultraviolet radiation^19–22,68 to above threshold or multiphoton ionization.^69–71 Our TDDFT approach explicitly includes only valance electrons, so we can only describe ionization processes related to valence electron dynamics.

Recently, density functional theory has been used in scattering calculations allowing the study of more complex molecules.^72–75 In the present paper, the dual propagation allows us to use wave packet propagation in large simulation boxes, and the scattering information can be extracted for any desired energy.

To solve the time-dependent Schrödinger-equation

we write the Hamiltonian as a sum of two parts,

where H₀ is the Volkov Hamiltonian describing the interaction of a free-electron with a laser field represented by a vector potential A(t) in velocity gauge,

and V(t) is the time-dependent potential. It is assumed that the time-dependent potential is only nonzero in a smaller region of the system.

The wave function will also be divided into two parts,

where ψ₀ is a solution of the TDSE of H₀,

Plugging ψ(r, t) into Eq. (1), we obtain an equation for ψ₁(r, t) as

where S is a source term defined as

Approaches like this have been often used in the solution of the time independent Schrödinger-equation (Lipmann–Schwinger method), and it has also been introduced to enforce outgoing wave boundary conditions in time-dependent systems.⁷⁶ 

The Volkov states, ϕ_k(r, t) (see  Appendix A), satisfy Eq. (5), so ψ₀ can be expanded as

Note that the c_k coefficients are time independent, so the time development of ψ₀(r, t) is described by the analytically known Volkov basis functions. These Volkov states can be used to solve the TDSE,^36,37 and the time propagation is very efficient using this basis, but one has to calculate the matrix elements of the potential, which becomes time consuming for very large systems.

Note that we can rewrite Eq. (6) in an alternative form as

with a source term defined by

In this case, the source term explicitly depends on ψ₁(r, t), which makes the solution more complicated (one has to use a predictor corrector step), but the Hamiltonian acting on ψ₁(r, t) is the Volkov Hamiltonian allowing “analytical” propagation. This alternative approach was not implemented numerically in the present work, but it might be useful in future applications.

In the calculation, the time propagation starts from an initial wave function ψ_initial(r), and at t = 0, we set

and

The initial state, ψ_initial(r), is the solution of the time independent problem with the vector potential A = 0. Equation (8) defines ψ₀(r, 0) in terms of the analytically known Volkov states, the source term S(r, t) can be calculated using Eq. (7), and Eq. (6) can be solved numerically with time propagation. Apart from possible numerical implementation inaccuracies, this approach has no approximations so far.

The source term is only nonzero in the region where the potential, V(r), is nonzero, so we can define a large enough region, C, that confines the interaction region. Region C will be used as our inner computational domain. When the time propagation of Eq. (6) starts, ψ₁(r, t) is zero, but due to the source term, it starts to evolve in time and space. After some time, ψ(r, t) approaches the boundary of C. If |ψ(r, t)|² is larger than a preset value, ϵ, at the boundary, we stop the propagation in the inner region. At this point, say, at t′ = t + ΔT, we expand the wave function into Volkov basis states as

and restart the propagation with

and

This step is repeated after every ΔT time interval. Note that Eq. (14) contains the whole wave function, describing the inner and outer regions, and nothing is lost or neglected.

The Volkov basis, ϕ_k(r, t), is defined in a large computational cell, which includes a sufficiently large part of the asymptotic region to describe an ionization process. The analytical time propagation of the wave function using this basis is not a computational burden, but it still requires a large memory storage for the c_k components. During the time propagation of ψ(r, t), one only needs to calculate the source term in region C and one only needs ψ₀(r, t) in C,

where $ϕkC(r,t)$ is a Volkov basis defined in region C. In this work, we have used Volkov states to time propagate ψ₀(r, t), but any other method would work.

The algorithm can be summarized as follows:

1. Initialize ψ₀(r, t) and $ψ0C(r,t)$⁠.

2. Propagate ψ₁(r, t) with Eq. (6) from t to t + ΔT.

3. Add ψ₀ and ψ₁,

1. Using Eq. (13), expand ψ into Volkov states ϕ_k(r, t + ΔT) and $ϕkC(r,t+ΔT)$ (defining ψ₀ and $ψ0C$⁠) and set ψ₁ = 0 and t = t + ΔT.

2. Go to step 2 and repeat the steps.

One can easily generate the expansion in Eq. (13) by fast Fourier transformation. This is the most expensive part of the calculation for large systems, but one only needs to do it in every ΔT time interval. The same expansion in terms of $ϕkC(r,t+ΔT)$ does not cause extra burden because region C is small.

Using the identity (“exponential time differencing”^45,77)

Eq. (6) can be rewritten as

This equation can be integrated from t to t + Δt,

From here, our goal is to numerically estimate the integral in Eq. (19) by approximating the source term.

We can approximate the integral in Eq. (19) assuming that source term S(r, t) is constant over the interval t to t = t + Δt,

which yields

where

This equation has an error of $O(Δt2)$⁠. Higher order approximations are also been tested.^45,77

To make the calculation practical, one has to represent the exponential operator. The simplest choice is a first order expansion,

and then

Alternatively, we can use Eq. (9) as a starting point for exponential time differencing. Repeating the steps above, using the first order expansion for the exponential, one can get a very simple equation for the time propagation with the $Ŝ$ source term,

because H_V is diagonal in momentum space and one can easily do the $e−iH0Δtψ1(r,t)$ part. The ψ₁ dependence of $S̃$ can be addressed with a predictor corrector step.

To increase the accuracy, one has to go beyond the first order approximation of the exponential. As we will implement the calculations on a real space grid, the simplest approach is a Taylor series expansion,

Using this expansion, Eq. (21) becomes

This propagation scheme is very similar to the conventional approaches and handles the wave function and the source term on the same footing. Note that the time step should be sufficiently small to preserve unitarity, but the approach proved to be stable and reliable in previous calculations,^47,78,79 and the small time step is also required by the timescale of dynamics of the system. This is the representation used in the calculations presented below.

We have tested the dual propagation on two different problems: (i) ionization of molecules and (ii) electron scattering on molecules. First, we show the results for one-dimensional tests and then for three-dimensional examples. We list all parameters (lattice spacing, box sizes, and time steps) for each example. These parameters are given in atomic units unless otherwise noted. The parameters are chosen after convergence tests, and the presented results are converged with respect to the time step and grid spacing.

In this section, we will use the Hamiltonian

where

is a soft Coulomb potential and q is a charge and W(x) is an absorbing potential to avoid spurious reflections from the boundary.⁸⁰ We adopt the following form for the complex absorbing potential (CAP):⁵⁴ 

Here, Δx = x₂ − x₁, x₁ is the start and x₂ is the end of the absorbing region, c = 2.62 is a numerical constant, and m is the electron’s mass. Function f is defined as

This form has been used in many calculations, and one of the advantages of this potential is that it has only one free parameter, the width of the potential. By increasing the width, the reflection can be decreased systematically.

First, we will consider a wave packet scattering on a potential with q = 1. If the momentum of the wave packet is high, then one needs very fine grid spacing due to the oscillatory nature of the wave function. To avoid that, we will separate the oscillatory part,⁸¹ 

with

and momentum of the wave packet p₀ = ℏk₀. This factorization is equivalent with setting A = p₀ in Eq. (28) and propagating ϕ_g(x, 0).

The initial wave packet is centered at x₀ = −10 a.u. with width a = 1 a.u. It is time-propagated up to T = 40 a.u. using the time step Δt = 0.005 a.u. and ΔT = 20Δt. ΔT has to be smaller than the time needed for ψ₁(r, t) to reach the boundary of the inner region. In the calculations, we monitor the density at the boundary, and |ψ₁(r, t)|² > 10⁻¹⁰ sets the upper limit on ΔT. After a short test, a suitable ΔT can be defined. In principle, one can also use variable ΔT adjusted by the density requirement at the boundary. Calculations were also repeated with different ΔT to test convergence.

The inner region, C, is from −20 to 20 a.u., and the larger region is from −250 to 250 a.u. with grid spacing 0.5 a.u. and a 50 a.u. wide CAP placed at both ends. Once the evolution of the wave packet is complete, one can transform the wave function from the time domain into energy space,

and the scattering information can be extracted at any desired energy. In Fig. 1, we show the results of the dual time propagation for two different momenta. We have also propagated the wave function in the full computational cell (from −250 to 250 a.u.) with the same time step, grid spacing, and CAP. We refer to this calculation as “full propagation.” The results of the dual propagation and the full propagation are in excellent agreement. We have also added the real part of the wave function calculated by the dual propagation to Fig. 1 to illustrate the oscillation of the scattering function.

The next example is the ionization of the one-dimensional hydrogen atom. In this case, we use q = −1 with two different laser pulses. The first one is defined as

with E₀ = 0.2, ω = 0.148, and τ = 50. This corresponds to a quiver distance of x_p = 2A₀/ωc = 18.3 a.u. and ponderomotive energy of $Up=A02/4c2$ = 0.37 a.u. Figure 2 shows the shape of this laser pulse. For the second pulse, we adapt the laser of Ref. 58. It is a linearly polarized laser pulse with λ = 532 nm (ω = 0.0856 a.u.) and peak intensity I = 1.38 × 10¹³ W/cm². The laser shape (see Fig. 2) is defined by

where f(t) is a trapezoidal envelope function of 14 optical cycles with two-cycle linear ramps, constant for ten cycles, and A₀ = 31.7 a.u., x_p = 5.4 a.u., U_p = 0.0133 a.u.

The results of these calculations are shown in Fig. 2. Figures 2(a)–2(f) show the time-dependent expectation value,

the probability amplitude |ψ(x, t)|², and the photoelectron spectrum (PES). The PES is defined as

where $ϕ̃(p)$ is the Fourier transform of (1 − M(x))ψ(x, t′) at time t′ after the duration of the laser pulse. Here, M(x) is the mask function (see  Appendix B) to cut off the inner part of the wave function and keep only the asymptotic part that contributes to the photoelectron emission. Different approaches to calculate the PES have been developed, compared, and discussed in Refs. 58, 60, and 64. The inner box is in the region of [−80, 80], and the large box is defined as [−600, 600], with Δx = 0.4 grid spacing and 100 wide CAP region (all in a.u.). The time step is Δt = 0.01 a.u. and ΔT = 25Δt. The first laser was propagated up to T = 300 a.u., the second up to T = 1200 a.u. We also compare the results to mask function calculations. The mask function starts at the boundary of the inner region to make the dual and mask method comparable, and the width of the mask function is 10 a.u. The same mask function is used as in Ref. 58 (see  Appendix B).

The first laser is a short strong few-cycle pulse, and the energies calculated by the dual and the full propagation are in perfect agreement. We also show the energy calculated in the inner region [using ψ₁(x, t)]. This energy differs from the energy calculated using the full wave function because of the contribution from the outer region. The second laser is a 14 cycle pulse with a smaller amplitude than that of the first laser. The energy change is smaller but very oscillatory, following the laser cycles. The agreement of the full and dual calculation is perfect.

Figures 2(c) and 2(d) compare the square of the wave functions after the laser pulse. The three approaches, the full, dual, and mask methods, are in perfect agreement. The slight disagreement close to the boundary is due to the difference in the application of the CAP in the mask and dual method and a slight scattering from the mask region in the case of the mask method. Figures 2(c) and 2(d) show that the shorter stronger pulse causes more ionization: a larger portion of the probability distribution is outside of the inner box, and the initial ground state wave function (it is a sharply peaked function at the origin) is barely noticeable after the pulse. The longer pulse leads the less ionization (the initial wave function is clearly visible) but steers the electron distribution farther away from the center.

The PES calculated by the two lasers is shown in Figs. 2(e) and 2(f). The three approaches are again in excellent agreement. The longer laser with more optical cycles seems to introduce more oscillations, but the general trend of decreasing intensity with increasing energy is similar.

These laser parameters are frequently used in test calculations,^58,59 but they are relatively weak. The dual propagation approach remains accurate for larger laser amplitudes (see Fig. 3), while the mask method starts to deviate from the full solution due to reflections from the boundary region.

As the mask region and region C is the same and the outer wave functions are updated the after ΔT = 25Δt time, the computational cost of the dual and the mask method is the same. The full propagation is updated at every Δt time steps, so the computational cost of the full propagation is about 25 times more than that of the dual propagation.

In this section, we show the application of the approach for electron scattering on benzene and uracil molecules and the calculation of the photoelectron angular distribution for an H atom and for CO and methyloxirane molecules. In the three-dimensional TDDFT calculations (including many electron orbitals), comparison to the full propagation is not possible, but we have carefully checked the convergence with respect to the time steps and grid spacings to test the approach. In the case of the H atom, there is only one orbital, and in this case, we will compare the dual propagation and the mask method.

In these calculations, we use the time-dependent density functional theory Hamiltonian to describe the electron dynamics. The TDDFT³⁸ is one of the most promising tools used to study the interaction of laser pulses with atoms and molecules.^{78,79,82–85} The TDDFT has been implemented using different basis sets, e.g., real space,^47,86 Gaussians,⁸⁷ or plane waves.^88,89 The TDDFT has been proven to produce accurate descriptions of the total and individual ionization yields for Ne and Ar atoms exposed to strong laser pulses⁸² and has also helped to explain the enhanced ionization in molecules as well as the energetics and dynamics⁹⁰ of laser-assisted field evaporation.⁹¹ Above threshold ionization and high harmonic generation have also been a subject of TDDFT studies.^{59,79,92–94}

In TDDFT, the electron dynamics is described by the time-dependent Kohn–Sham equation (TDKS),

where ψ_k represent fictitious non-interacting single particle orbitals, which yield the same density as the true electron wave function.

The Kohn–Sham potential, V_KS, may be decomposed as

Here, ρ is the electron density, which is defined by a sum over all occupied orbitals,

where the coefficient N_k accounts for the number of electrons in each orbital. V_H is the Hartree potential, defined by

which accounts for the electrostatic Coulomb interactions between electrons. V_XC is the exchange–correlation potential, and V_ion is the external potential due to the ions. The potential of the ions are represented by employing norm-conserving pseudopotentials of the form given by Troullier and Martins,⁹⁵ and the local-density approximation (LDA) approximation is used for the exchange–correlation potential.⁹⁶ 

For the ground state calculations, the molecular geometries were optimized by DFT calculations using the V_KS potential. In the calculations, the grid spacings are chosen to accurately represent the pseudopotential part of the V_KS potential. In the time-dependent calculations, the nuclear coordinates were kept fixed.

In the three-dimensional case, the PES is defined as

where $ϕ̃i(p)$ is the Fourier transform of $1−M(r)ψ(r,t)$ after the laser pulse. The photoelectron angular distribution is given by

Another experimentally useful quantity is the photoelectron circular dichroism (PECD) characterizing the asymmetry in the photoelectron angular distribution,

where θ is the polar emission angle of the electron with respect to the light propagation and P₊ and P₋ are the total PAD for the right and left polarized lasers, respectively. This quantity is mainly used in experiments measuring PAD of randomly oriented samples, but its dependence on molecular orientations has also been studied.³¹ 

As a first example, we study electron scattering on benzene and uracil molecules. The three-dimensional wave packet is defined as

with

The wave packet (a = 1 a.u.) is centered 19 a.u. away from the plane of the molecule, only the x component of its momentum is nonzero, and the molecule lies in the yz plane. In this case, only the wave function of the scattering electron is dual propagated, and the electrons of the molecule are time-propagated in the conventional way (Taylor propagator). The inner box has 61 × 61 × 61 points, and the large box has 305 × 305 × 305 points with 0.57 a.u. grid spacing. The time steps are Δt = 0.02 a.u. and ΔT = 10Δt. Figure 4 shows the scattered wave packets. Figures 4(a)–4(d) illustrate the interference created by the electron density distribution of the molecule. The difference between the uracil and the symmetric benzene clearly appears. Scattering of wave packets created by strong-field tunnel ionization has recently become an experimental tool,^25,27 and the dual propagation may help to simulate such experiments.

The photoelectron angular distribution of the H atom has been calculated by several research groups.^18,58,60,97 In these calculations, a laser pulse similar to

was used. In our calculation, we set the wavelength to λ = 800 nm (ω = 0.057 a.u.), and the intensity is I = 2.5 × 10¹⁴ W/cm², with N_c = 5 and A₀ = 91.3 a.u. In this example, the inner box is 101 × 101 × 101, the large box is 505 × 505 × 505, and the grid spacing is 0.4 a.u. The time steps are Δt = 0.02 a.u. and ΔT = 20Δt (x_p = 23.4 a.u., U_p = 0.11 a.u.).

The PAD is shown in Fig. 5. The radial distance from the origin corresponds to the electron energy, and the angle shows the direction of emission with respect to the laser polarization. The ring structure, corresponding to the above threshold ionization peak intensities, is very similar to the results of other calculations.^18,58,60,97 Using the same inner and outer regions, we have compared PAD calculated by the dual propagation (P_dual) and the mask method (P_mask). The difference |(P_dual − P_mask)/P_dual| is less than 10⁻⁴ for all energies and angles shown in Fig. 5.

Finally, we present a study of PES of molecules in circularly polarized light. In this case, all electronic orbitals are time-propagated. The left (+) and right (−) circularly polarized laser pulses have the form

and A(t) = −dE(t)/dt. The laser parameters are ω = 0.057 a.u., E₀ = 0.0836 a.u., t_c = 413.4 a.u., and τ = t_c/2. The inner box is 81 × 81 × 81, the large box is 405 × 405 × 405, and the grid spacing is 0.5 a.u. The time steps are Δt = 0.02 a.u. and ΔT = 20Δt.

Figure 6 compares the PES for a left and right circularly polarized pulse. Figure 6 also shows ∑ρ_i(p), the Fourier transform of the density, for comparison. The Fourier transform of the orbital density is defined as

where $ψ̃i(p)$ is the Fourier transform of orbital ψ_i(r, t) after the laser pulse. Figures 6(a) and 6(c) show the difference of the left steered and right steered electron density by the circularly polarized laser. As shown in Fig. 6, the left and the right polarized PAD are chiral images of each other in the present geometry.

The results for the C₃H₆O molecule are shown in Fig. 7. This calculation is a computationally expensive important test of the approach. In this case, 11 orbitals are propagated, and each orbital is represented on 405³ grid points in the outer box. In the inner box, each orbital is propagated using Eq. (6), and then the total electron density and potentials are calculated. Despite the size and the more complex coupled equations, the accuracy of the approach is maintained. We have checked different grid spacings and time steps, and the energy convergence and the final PAD remained the same. The test calculations have also shown that the present box side is suitable for the calculation. Much smaller boxes lead to very different PADs.

The momentum space densities (see Fig. 7) nicely show the effect of the left and right polarized lights in this case as well. Due to the more complex molecular structure, the PES for C₃H₆O is very different in the left and right polarization cases. Here, we have a very complex picture due to the orbitals with different binding energy and symmetry. Figure 8 shows the PECD as the function of the polar emission angle for the orientation shown in Fig. 7. The PECD asymmetry is about 2%, which is somewhat smaller than the 4% value measured for randomly oriented molecules.⁷¹ The calculations have to be repeated and averaged using many molecular orientations to get closer to the experimentally measured PECD.

To get further insights into these results, we picture the ground state orbitals in real and momentum space in Fig. 9. The orbitals have different symmetries in real space and in momentum space. The right and left polarized lasers ionize the orbitals to a different extent depending on the symmetry and shape and depth of the time-dependent potential induced by the laser field. Orbitals 7 and 9 ionized most in the right polarized laser, orbitals 10 and 11 ionized most in the left polarized laser (about 2% ionization in the present field). The ionization of orbitals 8 and 12 is about three times smaller, and the ionization of the lower orbitals is negligible. Figure 10 shows the momentum space density and PES for orbitals 7, 9, 10, and 11 after the duration of the laser. Figure 10 shows high momentum components of the orbital densities in momentum space. Comparing the momentum space densities before (Fig. 9) and after (Fig. 10) the laser excitation, one can see that the excited momentum distribution shows little resemblance to the ground state symmetry, and the right and left circularly polarized laser field strongly perturbed the original symmetry. Note that we had to use a different color scheme to amplify the differences, so the low momentum ground state components do not appear in Fig. 10.

Figure 10 shows that the major contribution to the PES of methyloxirane (Fig. 7) comes from orbitals 10 and 11 for the left and from orbitals 7 and 9 for the right polarized case. It is somewhat unexpected that the highest occupied molecular orbital (HOMO) is less ionized than the lower ones, but we have found similar cases previously.⁸⁴ In Ref. 84, we have studied the alignment dependent ionization of acetylene and ethylene by strong laser pulses, and we have found that the inner orbitals can be ionized to a larger extent than the HOMO. The ionization of the inner orbitals is primarily due to their extended weakly bound density tails.

We have implemented and tested a dual time propagation method that allows the simulation of ionization and electron scattering in large computational cells. The approach employs two wave functions, one in the inner region and the second one that is defined in a much larger space. The two wave functions are connected by a Schrödinger equation with a source term. Various one-dimensional ionization and scattering examples were used to test the accuracy of the approach by comparing a full solution to the dual propagation.

We have also presented three-dimensional examples. The electron scattering on molecules using wave packets allows the extraction of scattering information for a range of energies simultaneously. In this case, the wave packet has to be propagated in a large space starting from a position where the wave packet and the molecule do not interact and ending in a final position where the interaction is zero. During the propagation, the wave packet spreads in each direction requiring a large simulation box. The wave packet’s momentum is not constrained to be low, and for high momentum, we factorize out the oscillatory part that would require prohibitively small grid spacing. In principle, the approach can also be used for time-dependent scattering, e.g., when the molecule is subject to laser excitation and the wave packet is a probe, or the molecule’s excitation is caused by the wave packet itself.

The second set of examples illustrates the application for ionization. Ionization calculations obviously require large simulation regions, and the study of excitations by circularly polarized lasers further increases the size of the computational cells. We have studied the CO and methyloxirane PES in circularly polarized light. Photoelectron circular dichroism has been experimentally studied in methyloxirane.⁷¹ Photoelectron circular dichroism is normally studied in the gas phase, so the alignment of the molecule and the laser is random. In this work, we have calculated the PES for a given alignment, so direct comparison to experiment is limited, but measurements for aligned molecules will be possible in the near future.

The dual propagation greatly reduces the computational cost of the time-dependent calculations. The computational cost of the calculation in the inner box is the same as a conventional TDDFT calculation. The propagation in the large box requires a pair of fast Fourier transformation for each orbital. This Fourier transformation needs significant memory and computational time, but it has to be repeated only at about every 20 time steps. One can possibly make this step more efficient by using a nonuniform Fourier grid.

This work was supported by the National Science Foundation (NSF) under Grant No. IRES 1826917.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The Volkov Hamiltonian describing the interaction of a free-electron with a laser field represented in velocity gauge is given by

In this case, the TDSE admits a solution as Volkov states,

where

with

Within the mask method formalism,⁵⁸ the system is divided into two regions. The partitioning ψi(r) = ψ_A,i(r) + ψ_B,i(r) is enforced by

where M(r) is a mask function that smoothly connects the two regions. In region (A), the time evolution is performed with the full Hamiltonian. In region (B), the electron–nuclear and electron–electron interactions are neglected, and thus, the wave function may be projected into momentum space via a fast Fourier transform (FFT) and propagated analytically using the velocity gauge Volkov expansion. Another advantage of the mask method is that the mask function does not need to be enforced at every time step, and this mask time step is usually defined as h = nΔt. The time propagation for the full wave function proceeds as follows:

1. Propagate the real-space component ψ_A(r, t) by applying the time evolution operator (or another method) to time t + h.

2. Propagate the momentum space component $ξB(k,t)=F{ψB(r,t)}$ by applying the Volkov propagator to time t + h.

3. Reverse Fourier transform ξ_B(k, t + h) into ψ_B(r, t + h) and convert to the length gauge.

4. Enforce the mask function to mix regions A and B.

5. Convert back to the velocity gauge and Fourier transform ψ_B(r, t + h) to momentum space ξ_B(k, t + h).

6. Repeat steps 1–5 until time t_final is reached.

In this paper, the mask function is defined to be

as used in Ref. 58. Because the real-space wave function in region B is known at times when the mask function is enforced, a position based mask can also be used to prevent artificial wrapping of the wave function into the opposing side of the simulation box.