Strong-field-induced single and double ionization dynamics from single and double excitations in a two-electron atom

Two-electron dynamics of an excited model atom interacting with moderately strong laser fields is analyzed in the time domain. We solve the time-dependent Schrödinger equation (TDSE) for two electrons confined to the same one-dimensional configuration space, accounting also for the electron-electron interaction. The computational method allows direct access to the time-dependent population of the relevant atomic states during and right after the interaction with a near-infrared (NIR) laser pulse. We compare the ionization dynamics of singly excited states and doubly excited states. We find that doubly-excited initial states exhibit enhanced double ionization yield, with non-trivial dynamics including contributions from direct and sequential processes, while the electrons leave the atom either back-to-back or in the same direction.


Introduction
The first step towards developing a complete picture of the intricate dynamics in multi-electron systems is to understand from bottom up the correlated two-electron atoms and molecules. The simplest atomic system of this kind, complex enough for its (classical or quantum) equations of motion not to be analytically solvable, is the helium atom. For this reason, helium has become a benchmark system to study electron-correlation effects in their great complexity, both theoretically and experimentally [1][2][3][4][5][6][7]. Ultimately, after gaining some knowledge about the system, one would like to influence the electron dynamics upon interaction with strong fields, which are an ideal tool to control the two-electron dynamics at the atomic level.
Typical experiments for studies of the electron dynamics in helium on the natural attosecond time scale group in two classes: all-optical techniques, such as transient absorption spectroscopy (TAS) [8] and methods to detect fragments after ionization. The latter include electron-detection experiments such as time-of-flight (TOF) [9] or velocity-map-imaging (VMI) setups [10] and full-coincidence experiments [11] carried out with reaction microscopes (ReMi) or cold target recoil ion momentum spectroscopy (COLTRIMS). The choice of the experimental technique is tightly linked to the nature of the studied electron process of interest. The all-optical TAS can access bound-state dynamics and couplings between states by identifying absorption lines and their changes upon near-infrared (NIR) light-matter interaction [6,[12][13][14][15]. However, this method gives only limited access to the NIR-driven electron dynamics during and after ionization. Electron-detection experiments, on the other hand, are used to capture the ionization dynamics through the direct detection and characterization of the outgoing electrons [11,[16][17][18]. For a complete characterization of the ionization process, full-coincidence experiments are the method of choice, where one is able to measure all possible ionization fragments, i.e., all electrons and ions, and quantify the yield and properties of the out-going particles in a kinematically complete experiment.
Theoretical investigations on helium which numerically solve the time-dependent Schrödinger equation (TDSE) either in its full dimensionality [3] or for models with reduced complexity [19][20][21], or specifically address the structure of the excitation spectrum [22], examine the effect of light-matter interaction on observables that are accessible in the previously mentioned experiments. The effects of strong NIR fields on the Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. electron dynamics, in this context, have been typically studied starting from the ground state influenced by a single NIR laser pulse [19,20] or a combination of an NIR and a time-delayed extreme ultraviolet (XUV) pulse, as in TAS [23][24][25]. One-and few-photon double ionization of helium has also been investigated under the influence of intense XUV fields, where couplings between the ground state and single and double excitations were observed [26,27].
Here, we concentrate on the electron dynamics in photo-excited states (i.e., accessible via dipole-allowed transitions from the ground state), considering both single and double excitations as the initial states, and compare the occurring ionization mechanisms under the influence of moderately strong NIR fields. We find both single and double ionization to play a role in strong-field ionization in both considered cases. Double ionization, however, is enhanced for doubly excited states (DESs) compared to singly excited states (SESs). We find both direct and sequential processes to contribute to double ionization during the interaction with the NIR field, which is revealed in the probability density of the doubly excited state in position space. The NIR intensity regime considered in this work matches typical intensities in previous experiments, where doubly excited resonances in helium vanish in absorption spectra [6,14] for~-I 10 W cm NIR 13 2 . For lower NIR intensities, on the order of -10 W cm 12 2 , it was found in these experiments that both singly and doubly excited resonances experience the same effect from the NIR pulse, in particular a light-induced phase shift resulting in an altered resonant absorption line shape [13]. While in the latter work laser-induced bound-bound transitions in helium have been compared across singly and doubly excited states, we now investigate how bound-free transitions, i.e., the electron ionization dynamics, differ for single and double excitations under the influence of a moderately strong NIR laser field.

Methods
To investigate the strong-field ionization dynamics from single and double excitations in helium we numerically solve the time-dependent Schrödinger equation for two electrons in a linearly polarized laser field, which allows to restrict the motion of each of the electrons to a one-dimensional (1D) discretized grid. We then separate the electronic and nuclear degrees of freedom, and do not further consider the motion of the nucleus. Since we solve the non-relativistic TDSE, we neglect the magnetic-field component, having an amplitude two orders of magnitude smaller than the electric field ( = = E B c 137; 0 0 / in atomic units). The described simulation approach is based on a non-perturbative two-electron model [28,29] and accounts for the full electron-electron correlation, as well as for linear and non-linear light-matter interaction under the influence of the NIR laser.
We use the non-relativistic Hamiltonian in the length gauge of the dipole approximation (in atomic units) Each of the Coulomb potentials is modified by a soft-core parameter, a for the electron-nucleus attraction and b for the electron-electron repulsion. This form of the softened Coulomb potential is known as the Rochester onedimensional potential [30,31]. All results presented in section 3 were obtained for = = a b 1 a.u. The NIR laser field is approximated by an oscillating cosine wave under a Gaussian envelope in the form E 0 is the maximum of the electric field strength, ω is the laser frequency, f is the carrier envelope phase of the pulse and t G is the Gaussian width of the pulse and is therefore associated with the duration of the laser pulse t P , defined as the full width at half maximum (FWHM) of E t .
2 | ( )| The connection between t G and t P evaluates as t t = 2 ln 2 .

P G
The time-dependent wave function of the model system in position space is denoted by y x x t , ; .  ( ) The absolute value squared of the complex number entries gives the electron probability density for each position.
The coordinates x 1 and x 2 span between -1023.5 a.u. and 1024 a.u. with D = x 0.5 a.u.
i To avoid the propagation of the wave function over the edge of the grid, and in order to simulate ionization effects, a cylindrically symmetric imaginary potential , V acting as an absorbing boundary, is employed to dampen the electron probability density at the outer regions of the grid. Figure 1(a) shows the grid with the absorbing boundary as shaded area. For more details on the numerical implementation of the position-space wave function and the potential V please refer to appendix A.
The described grid is partitioned into different sections, figure 1(a). This partitioning technique has been successfully implemented previously in several theoretical studies of the ionization from the ground state of small atoms and molecules [32][33][34][35][36][37][38]. Time-dependent changes in the electronic density in each of the grid segments, I to IX in figure 1(a), indicate dynamical effects caused by the interaction with an external laser field.
How the electronic density propagates in the different regions is especially influenced by the correlation between the electrons. The partitioning of the wave function is done in such a way that section V defines the size of the atom, meaning that at least 99% of the unperturbed two-electron probability density of the considered singly and doubly excited states is contained in this segment. Figures 1(b) and (c) show an atomic SES and DES, respectively. The boundary between this neutral atomic region and the neighbouring regions, i.e. the atomic radius, is at = d 64 a.u. atom Observing an increase of the electronic density in sections II, IV, VI and VIII (yellow sections) means that one of the electrons is predominantly located near the nucleus, whereas the other electron moves away from it. This is the process of single ionization (SI), where one of the electrons leaves its bound state, hence these four regions are termed SI sections. If the second electron also occupies regions away from the origin, captured in sections I, III, VII and IX (green and blue sections), both electrons leave the atom and one observes double ionization (DI), so these regions are called the DI sections. Furthermore, one can also distinguish between a strongly directional two-electron emission pattern [17,39,40] in sections III and VII, and a back-to-back emission pattern [16,40] in section I and IX. The goal of this work is to study the difference between ionization from a singly and a doubly excited state, when both are being exposed to a moderately strong NIR laser intensity. We therefore focus on the investigation of the NIR strong-field ionization dynamics from two selected singly and doubly excited states. The NIR pulse hereby is centered at photon energy w = =  1.74 eV 0.064 a.u., with duration of t = 5 fs P and intensity ranging from 0 to -14 TW cm 2 . Such a scenario could be experimentally realized after a preceding population of the excited states from the ground state by narrow-band XUV radiation, transferring electronic population from the ground state only to a specific excited state. In the simulation, we prepare the excited bound states with known energies to be excited eigenfunctions of the free Hamiltonian of the model system. For this we use the methods presented in [41] and in more detail in appendix B. Hereby, we deliberately exclude the initial XUV excitation step and therefore any laser interaction with the ground state of the model system. To compare the ionization behavior of the considered SES and DES we chose them according to their single-ionization potential I p , which for both states is around 0.1 a.u., more specifically 0.081 a.u. for the SES and 0.114 a.u. for the DES. A level scheme with the energy values for excited states of the field-free model system is provided in figure 1(d). The time evolution of the states is found as a solution of the TDSE, which is solved via a split-step algorithm method with second order accuracy [42].
3. Comparative analysis of the ionization dynamics from single and double excitations 3.1. Ionization probability In order to study ionization, we analyze the time-dependent state probability density y y = á ñ y P t t t ( ) ( )| ( ) . As mentioned in the previous section, only the central grid part, section V, contains the electronic bound states. The temporal evolution of y P t ( ) for this grid area is calculated and shown in figures 2(a.2) and (b.2) for the SES and the DES, respectively. The field strength E 0 of the NIR laser varies between 0 a.u. and 0.0225 a.u., corresponding to intensity variation in the range --0 14 TW cm 2 . The probability density of the unbound fraction is distributed over all other grid sections and is hence defined as the ionization probability =y P t P t 1 I ( ) ( ). We note that the so defined ionization probability depends on the exact choice of the box partition. 2) Ionization probability P I and state probability density y P t ( ) for a doubly excited state y . DES In both cases y P t ( ) is calculated only for region V in figure 1(a), where the electronic density in all other sections and the parts absorbed by the damping potential define the ionization probability P I for each state. All data is calculated by scanning E 0 linearly in steps of 0.0025 a.u. Note that the intensity scale is increasing quadratically with the field strength. Figures 2(a.1) and (b.1) depict P I right after the end of the laser pulse, at = t 500 a.u., as a function of NIR field strength for the SES and DES, respectively. Hereby, in the case of the DES, we correct for the autoionization rate of the states by subtracting the field-free single ionization yield, which remains below 0.4%. For both considered excited states, the ionization probability increases with increasing intensity and for~-I 10 W cm 13 2 approaches 50%, meaning that half of the two-electron probability density is no longer in the original state, i.e., ionization gets increasingly likely. The Keldysh parameter g = 2 , for both SES and DES is larger than 1 in the considered moderate NIR-intensity regime, which points towards a multi-photon ionization process. Using the condition for channel closing w +  m I U P P , we expect that at = E 0.015 a.u. 0 , which is marked by the horizontal dashed line in figures 2(b.1) and (b.2), the two-photon ionization channel for the doubly excited state is closed by the rising ponderomotive energy shift due to the electric field strength. The three-photon ionization channel then becomes the dominant one, as observed in experiments with noble gases in the strong-field limit [43,44]. Channel closing, caused by ponderomotive energy shifts is thus a candidate for explaining the slightly decreasing ionization probability of the DES for NIR intensity above -10 TW cm 2 .

Time evolution of the excited states
Considering that a significant amount of the electron density of both states spreads out of the neutral atom region for increasing laser intensity (figure 2), we now examine the wave-function dynamics in these outer ionization regions in more detail in For the singly excited state, the electrons occupy almost exclusively the regions where one of them is closely bound to the nucleus and the other one moves away from the neutral atom, i.e., the SI sections. A very small part of the electron probability density, on the level of the numerical precision, leaves the SI regions and propagates towards the DI regions. The evolution of the doubly excited state y , DES 2 | | on the other hand, looks different. Shortly before the maximum of the laser field, at time point (1) in figure 3(b), again the single ionization regions are most dominantly populated. For later times, however, a considerable amount of y DES 2 | | (note that the color scale for both states and for all times is kept same) enters the diagonal regions, where both electrons simultaneously appear away from the neutral atom region. In the back-to-back region I and IX, the electronic probability density moves perpendicularly out of the SI regions, which is manifested as rectangular structures in the DES time evolution. This is an indication of an ionization mechanism, where first one of the electrons is driven away from the bound state, i.e. ionized, and afterwards also the other electron moves to the DI region, also leaving the bound state. Horizontal or vertical flow of the wave function perpendicularly out of the closest SI region is therefore a signature of an independent, thus sequential double-ionization process. In the unidirectional regions III and VII, the electrons move diagonally either out of the central (bound) region or out of the nearest SI domains (indicated by black arrows in figure 3(b)) or also perpendicularly out of the SI regions, as in the back-to-back case. Hereby, the diagonal motion, where both electrons move simultaneously out of the neutral region, is an indication of a direct double ionization pathway. In many earlier studies, recollision-induced non-sequential ionization is found to play a significant role in strong-field ionization out of the atomic or molecular ground state [17][18][19]34], where the NIR intensity is on the order of -10 W cm 14 2 . In this study, however, we focus on the effects of moderately strong fields, reaching up to -10 W cm 13 2 . The maximum energy of a recolliding electron then amounts to = = U 3.17 0.076 a.u. 2.07 eV P , which is not sufficient to ionize the remainig bound electron, with binding energy of = 1.5664 a.u. 42.6 eV for the singly excited state and » 0.5 a.u. 13.6 eV for the doubly excited state. Recollision-induced double ionization is therefore not expected to play a dominant role for the here considered NIR intensity.

Quantification and directionality of single and double ionization
To quantify the population of the SI and DI regions we make use of the partition technique and determine the single and double ionization yield (SIY and DIY) for both considered states, through integrating over the probability density in the respective regions. As noted above, the so defined yields again depend on the exact choice of the box partition. According to our choice (see figure 1), the results are plotted in figure 4, where (a.1)  and (a.2) show the time evolution of the SIY and DIY, respectively, for the singly excited state, as a dependence on the maximum field strength of the NIR pulse. The same plots for the doubly excited state are shown in figures 4(b.1) and (b.2).
Due to the autoionizing nature of the doubly excited state, its initial SIY is higher than the initial SIY for the SES. In the field-free case the SIY of the DES decreases, again due to autoionization. For both states the increase in single ionization arises at earlier times compared to the double ionization, where the latter sets in around the overall maximum of the laser pulse, which is also confirmed by the evolution of the spatial wave function, shown in figure 3. The observation of the DIY rising after an initial increase of the SIY supports the above statement that a sequential process contributes to the observed double ionization. This is also in agreement with the results plotted in figure 3, where in both the back-to-back regions (I and IX) and the unidirectional regions (III and VII), perpendicular ionization structures are present, indicative of a sequential double ionization mechanism. The double ionization yield of the doubly excited state is increased by four orders of magnitude during the interaction with the NIR laser field, which is not surprising considering that the DES exhibits a lower double ionization threshold because the second electron is left in an excited state after single ionization and therefore it is more loosely bound. For the SES, we identify single ionization as the dominant ionization channel, as is evidenced by the small DIY for all electric field strengths.
We further analyze the contributions to the double ionization yield of the DES, distinguishing now between two double ionization scenariosback-to-back and unidirectional ionization. To compare the ionization yield for both cases, we integrate over regions I and IX for the back-to-back case and regions III and VII for the unidirectional double ionization, depicted in figures 5(a.1) and (b.1), respectively. For all considered NIR field strengths the back-to-back DIY is higher than the unidirectional yield. This is expected to be a consequence of the electron-electron interaction, where the unidirectional channel is suppressed through the Coulomb repulsion between the electrons. We track the motion of one electron in the double ionization regions by integrating the electron probability density over the distribution of the other electron. For the data shown in figure 5 we integrate over the x 2 -axis as indicated by arrows in all four DI regions. Figures 5(a.2)  Bursts of electronic probability density enter the double ionization regions, when the NIR field is the strongest. In both back-to-back and unidirectional regions the electrons are then driven by the oscillating laser field. There are a few different signatures between the back-to-back and the unidirectional cases. First, the electron probability density reaches the unidirectional region later than the back-to-back region, i.e., the back-to-back ionization starts before the unidirectional ionization. Second, the spread of the x 1 -projection in the back-to-back region is larger than in the unidirectional region. Furthermore, electrons escape the atom in the same direction predominantly around the envelope maximum of the electric field, as indicated by the lack of electron bursts for later times in figures 5(b.2) and (b.3) near the nucleus. Higher laser intensities are thus necessary for the electrons to overcome the additional Coulomb repulsion when escaping into the regions III and VII.

Conclusion
In this numerical study, we explored the role of the initial state in light-matter interaction by comparing the strong-field ionization dynamics of single-and two-electron excited atomic states in the limit of moderately strong fields. While single ionization is dominant for both single and double excitations, we observe characteristic differences in the double ionization behavior. After the onset of single ionization during an ultrashort NIR laser pulse, doubly excited states also exhibit significant double ionization. In our analysis we find both electrons escaping the nucleus either on the same side (unidirectional) or in the opposite directions (backto-back), where the latter pathways dominate the emission pattern. Furthermore, a sequential ionization process contributes to double ionization both for back-to-back as well as unidirectional emission, whereas a direct double ionization mechanism can be observed in addition for unidirectional emission. The observed enhanced double ionization yield only for initially doubly excited states could therefore be a key component to understanding the overall different ionization dynamics of single and double excitations. Excited states, in this sense, provide a new platform for double and multiple ionization studies to understand sequential and direct mechanisms beyond recollision. We note that experimental measurements of such dynamics are feasible by employing coincidence detection of both electrons. To this end, using the selective excitation of different singly and doubly excited states by means of XUV excitation, the subsequent strong-field-ionization pathways can be studied with an NIR laser pulse following an XUV pulse. With the availability of tunable narrowband XUV sources at free-electron lasers [45,46], a significant population of different excited states is possible [47], enabling the experimental measurement of strong-field ionization dynamics as proposed here, while being sensitive to the role of the electron correlation within the initial state. The absolute value of the boundary is shown in figure 1(a) as a shadowed area.

Appendix B. Preparation of the excited states
Before we prepare the excited states of the model atom, we first consider the properties they should have, with respect to symmetry of the wave function. To account for Pauli's exclusion principle, we are hereby only interested in wave functions that are antisymmetric with respect to electron exchange and therefore in combinations of spatially symmetric and spin antisymmetric wave functions or vice versa. This way we intrinsically include spin correlation in our model. Since transitions from singlet (spin-antisymmetric) to triplet (spin-symmetric) states are dipole forbidden, we remain in the singlet manifold at all times, i.e., all states' spatial wave functions are symmetric with respect to electron exchange. For this reason, we first find the eigenenergies of all spatially symmetric eigenfunctions of the model atom. We can use two different methods to find the eigenenergies of the states of the model atom: (1) by first finding the ground state of the system and its energy by means of negative imaginary time propagation [48] and then exciting transitions from the ground state to excited states with an XUV pulse. From the energy of the resonance transitions in the absorption spectrum we then obtain the energy of the excited states; (2) by a spectral method where we use a spatially symmetric random wave function, which we propagate in time and at every time step we evaluate the inner product of the wave function at time t with the initial random wave function, which results in a set of time-dependent projection coefficients a t ( ). The Fourier transform of a t ( ) includes all eigenenergies of the model atom. The results of both methods agree with each other and we therefore obtain the eigenenergies of the spatially symmetric excited states of the model atom, shown in figure 1(d). Precisely knowing the energy of an excited state, we can prepare it by multiplying a random symmetrized grid representation of the overall 2D spatial wave function at each time of its free time evolution with the factor w e i t m , where w m is the energy of the state we wish to extract. Writing y x t , ( )as a combination of all eigenstates of the system, we then obtain Besides the state y m š we wish to extract, all other states acquire an additional phase shift, so on average the contribution of those states can be neglected for a sufficiently long time propagation. On the other hand, at each time, the contribution of the state with energy w m is effectively enhanced, so that we end up with the desired excited state of the model atom. To crosscheck the accuracy of this method, singly excited states were also prepared using the alternative approach of negative imaginary time propagation when successively projecting out all bound states with smaller energy, as suggested in [48]. We find the two methods to show identical results for the calculation of electronic probability, within the numerical precision of -10 6 .