A general model and toolkit for the ionization of three or more electrons in strongly driven molecules using an effective Coulomb potential for the interaction between bound electrons

We formulate a general three-dimensional semiclassical model for the study of correlated multielectron escape during fragmentation of molecules driven by intense infrared laser pulses, while fully accounting for the magnetic field of the laser pulse. We do so in the context of triple ionization of strongly driven HeH$_{2}^{+}$. Our model fully accounts for the singularity in the Coulomb potentials of a recolliding electron with the core and a bound electron with the core as well as for the interaction of a recolliding with a bound electron. To avoid artificial autoionization, our model employs effective potentials to treat the interaction between bound electrons. We focus on triple and double ionization as well as frustrated triple and frustrated double ionization. In these processes, we identify and explain the main features of the sum of the kinetic energies of the final ion fragments. We find that frustrated double ionization is a major ionization process, and we identify the different channels and hence different final fragments that are obtained through frustrated double ionization. Also, we discuss the differences between frustrated double and triple ionization.


I. INTRODUCTION
Multielectron ionization and the formation of highly excited Rydberg states are fundamental processes in molecules driven by intense and infrared laser pulses.Rydberg states have a wide range of applications, such as, acceleration of neutral particles [1], spectral features of photoelectrons [2], and formation of molecules via longrange interactions [3].Electron-nuclear correlated multiphoton resonance excitation was shown to be the mechanism responsible for the formation of Rydberg states in weakly-driven H 2 [4].This process was shown to merge with frustrated double ionization for H 2 driven by intense infrared laser fields-strongly driven [5,6].In frustrated ionization an electron first tunnel ionizes in the driving laser field.Then, due to the electric field, this electron is recaptured by the parent ion in a Rydberg state [7].
Most studies on strongly driven molecules address double and frustrated double ionization of two-electron molecules [5,6,[8][9][10][11][12][13][14][15][16][17][18].However, there are scarcely any theoretical studies on three-electron ionization and formation of Rydberg states during the fragmentation of strongly driven molecules.The reason is that accounting for both multielectron and nuclear motion is a formidable task.This is corroborated by the few theoretical studies for three-electron ionization of strongly driven atoms, a simpler but still highly challenging task.Mostly formulated in the dipole approximation, these studies on atoms employ lower dimensionality classical [19,20] and quantum mechanical [21,22] models to reduce the complexity and computational resources required.However, lower dimensionality results in a non accurate description of electron-electron interaction during triple ionization.Currently, only classical or semiclassical 3D models of triple ionization of atoms are available [20,[23][24][25][26].
In Refs.[27,28], we argued that the main disadvantage of available classical and quantum models of triple ionization in atoms is their softening of the Coulomb interaction of each electron with the core.In quantum mechanical models, this interaction is softened to obtain a computationally tractable problem.However, in classical and semiclassical models, softening the Coulomb singularity is fundamental and relates to unphysical autoionization.Classically there is no lower energy bound.Hence, when a bound electron undergoes a close encounter with the core, the Coulomb singularity allows this electron to acquire a very negative energy.Through the Coulomb interaction, this energy can be shared by another bound electron potentially leading to its artificial escape.
To avoid artificial autoionization, most classical and semiclassical models of triple ionization of strongly driven atoms soften the Coulomb potential [20,23,24] or add Heisenberg potentials [29] (effective softening) to mimic the Heisenberg uncertainty principle and prevent each electron from a close encounter with the core [25,26].However, softening the Coulomb potential does not allow for an accurate description of electron scattering from the core [30,31].This is due to the exponential decrease of the ratio of the scattering amplitude for the soft-core potential over the one for the Coulomb potential with increasing momentum transfer [30,31].For recollisions [32], this implies that soft potentials are quite inaccurate for high energy recolliding electrons that backscatter.Indeed, we have shown [27,28] that the ionization spectra obtained with the Heseinberg model differ from experimental ones obtained for driven Ne and Ar [33][34][35][36][37][38][39].
Here, we formulate a general 3D semiclassical model of multielectron ionization and fragmentation of strongly driven molecules, while fully accounting for the magnetic field of the laser pulse.The main premise in our model is that two interactions are most important during a recollision or return of an electron to the core and hence are treated exactly.We account for the singularity in the Coulomb potential between each electron, bound or quasifree, and the core.Quasifree refers to a recolliding electron or an electron escaping to the continuum.Also, we treat exactly the Coulomb potential between each pair of a quasifree and a bound electron and hence the transfer of energy from a quasifree to a bound electron.However, accounting for the singularity in the electron-core interaction, implies we need to avoid unphysical autoionization that can take place through energy transfer between bound electrons.We do so by approximating the energy transfer from a bound to a bound electron via the use of effective Coulomb potentials.We refer to this model that accurately accounts for all Coulomb interactions and only employs approximate effective Coulomb potentials to describe the bound-bound electron interaction as ECBB model.This is a sophisticated model that identifies during time propagation whether an electron is quasifree or bound.That is, in the ECBB model, we decide on the fly if the full or effective Coulomb potential describes the interaction between a pair of electrons.We do so by using a set of simple criteria detailed below.
The ECBB model for strongly driven molecules is a generalization of the ECBB model we have previously developed for strongly driven atoms [27,28].For atoms, we have shown that the ECBB model results in triple ionization spectra of strongly driven Ne that are in excellent agreement with experiment [28].The main difficulty in generalizing the ECBB model to molecules is the formulation of the effective Coulomb potential between a pair of bound electrons in the presence of many nuclei.In atoms, the effective Coulomb potential is obtained by assuming that a bound electron has a charge distribution centered around the atomic core.In molecules, we assume that an electron creates an effective potential that is the sum of the effective Coulomb potentials with respect to each nucleus separately, i.e., each effective Coulomb potential is the same as for an atom.However, each atomic effective Coulomb potential is weighted by an approximate expression for the probability to find the electron that creates this potential at a given position around this nucleus.
Using the ECBB model, we address triple and double ionization in the triatomic molecule HeH + 2 when driven by a linearly polarized laser pulse.We also address the formation of Rydberg states via frustrated ionization.Employing the HeH + 2 molecule allows us to directly compare the results for triple, double and frustrated triple ionization obtained with the ECBB model in this study and with a predecessor of the ECBB model obtained in previous work [40].This previous model determined on the fly whether an electron is quasifree or bound, as the ECBB model does, but this was done using a set of criteria less sophisticated than the ones employed in the ECBB model for atoms [27,28] and molecules in the current work.Also, in this previously employed model the interaction of a pair of bound electrons was set equal to zero.As a result, in Ref. [40], we could only address frustrated triple ionization where two electrons ionize and one electron remains bound in a Rydberg state.Here, since we account via effective potentials for the interaction between bound electrons, we can also address frustrated double ionization.In this process, one electron ionizes, the other remains bound in the ground state of one of the He or H fragments and another electron remains bound in a Rydberg state of one of the fragments.
Here, we address both triple and double ionization (TI and DI) as well as frustrated triple and double ionization (FTI and FDI) and obtain the sum of the kinetic energies of the atomic fragments, referred to as kinetic energy release.We find that the atomic fragments are moving faster when fragmentation of HeH + 2 is described with the ECBB model compared to its predecessor [27,28].This is consistent with our finding that the electron-electron escape to the continuum is more correlated when the process is described with the ECBB model, since it allows for energy transfer between bound electrons.Higher correlation results in faster electron escape leading to Coulomb explosion of the nuclei at smaller distances, eventually resulting in higher kinetic energies.Also, we find that, in three-electron molecules, FTI and FDI proceed via two pathways, first identified in FDI of the strongly driven two-electron molecule H 2 [6].One electron ionizes early on (first step), while the remaining bound electron does so later in time (second step).If the second (first) ionization step is frustrated, we label the FTI and FDI pathway as A and B, respectively.

II. METHOD
In what follows, we describe in detail the formulation of the ECBB model for strongly driven molecules.The ECBB model resolves unphysical autoionization in 3D semiclassical models that fully account for the Coulomb singularity.Also, we formulate the ECBB model in the nondipole approximation fully accounting for the magnetic field component of the laser field.Finally, both electrons and the cores are propagated in time.

A. Definition of the effective charge and potential
In what follows, the cores are assigned indices in the interval [1, N c ], where N c is the number of cores.The remaining indices in the interval [N c + 1, N ], with N the number of particles, are assigned to the electrons.For each electron j, we define an effective charge ζ j,n (t) associated with each core n as follows [41] where Q n is the charge of the core n, E n 1s is the groundstate energy of a hydrogenic atom with core charge Q n , i.e., E n 1s = −Q 2 n /2, and E j (t) is the energy of electron j given by where A(r j , t) is the vector potential and E(r j , t) = − ∂A(rj ,t) ∂t is the electric field.The effective potential that an electron j experiences at a distance |r n − r j | from the core n due to the charge of electron i is given by [41] V To generalize this effective Coulomb potential to molecules, in Eq. ( 2), we assume that an electron i screens each core n seperately as if it was an atom, with probability |C i,n | 2 given by with Hence, we approximate the wavefunction of each bound electron i with a 1s hydrogenic wavefunction around each core.In addition, we are distributing the charge of a bound electron i amongst the different cores according to the probability density of electron i with respect to each core.Hence, the total potential that a bound electron j experiences due to a bound electron i is given by

B. Hamiltonian of the system
The Hamiltonian of the N -body molecular system, comprised of N c cores and N -N c electrons in the nondipole approximation is given by where Q i is the charge, m i is the mass, r i is the position vector and pi is the canonical momentum vector of particle i.The mechanical momentum p i is given by The potential V i,j is given as a sum of effective potentials as follows The functions c i,j (t) determine whether the full Coulomb interaction or the effective V eff (ζ i,n (t), |r n − r j |) and V eff (ζ j,n (t), |r n − r i |) potential interactions are on or off for any pair of electrons i and j during the time propagation.Specifically, the limiting values of c i,j (t) are zero and one.The value zero corresponds to the full Coulomb potential being turned on while the effective Coulomb potentials are off.This occurs for a pair of electrons i and j where either electron i or j is quasifree.The value one corresponds to the effective Coulomb potentials V eff (ζ i,n (t), |r n − r j |) and V eff (ζ j,n (t), |r n − r i |) being turned on while the full Coulomb potential is off.This occurs when both electrons i and j are bound.For simplicity, we choose c i,j (t) to change linearly with time between the limiting values zero and one [27,28].Hence, c i,j (t) is defined as follows where c(t) = β(t − t i,j s ) + c 0 , and c 0 is the value of c i,j (t) just before a switch at time t i,j s .A switch at time t i,j s occurs if the interaction between electrons i, j changes from full Coulomb to effective Coulomb potential or vice versa.Every time during propagation that such a switch takes place, we check whether for each pair of electrons the full Coulomb force should be switched on and hence the effective potential switched off or the full Coulomb force should be switched off and the effective potential switched on.The former occurs if at time t i,j s one of the two electrons in a pair of bound electrons changes to being quasifree while the latter occurs if in a pair of a quasifree electron and a bound electron the quasifree electron becomes bound.At the start of the propagation at time t 0 , t i,j s is equal to t 0 and c 0 is one for pairs of electrons that are bound and zero otherwise.To allow the effective Coulomb potential to be switched on or off in a smooth way, we choose β equal to ±0.1; plus corresponds to a switch on and minus to a switch off of the effective Coulomb potential.

C. Global regularisation
We perform a global regularisation [42] to avoid any numerical issues arising from the Coulomb singularities.We previously used this regularisation scheme, for strongly driven H 2 , to study double and frustrated double ionization within the dipole approximation [43] as well as nondipole effects in nonsequential double ionization [44].Also, we have used this regularisation to study triple ionization in the nondipole approximation in strongly driven Ar [27] and Ne [28].In this scheme, our new coordinates involve the relative position between two particles i and j q ij = r i − r j (11) and their conjugate momenta where The inverse transformation is given by and where Next, we define a fictitious particle k for each pair of particles i, j as follows with j > i and the total number of fictitious particles being equal to K = N (N − 1)/2.In addition, we define the parameters α ik and β ik , as Given the above, Eqs. ( 14) and ( 15) take the following simplified form and

D. Derivation of the time derivative of the effective charges
The Hamiltonian in Eq. ( 7) depends not only on positions, momenta, and time but also on the effective charges.Moreover, the Hamiltonian depends on time through the vector potential as well as through the effective charges that are time dependent.Since the effective charge of electron j is proportional to the energy E j (t) [see Eq. ( 2)], it follows that we must obtain the derivative with respect to time of E j (t).This is necessary at any time during propagation if at least two electrons are bound.Following the same procedure as in Ref. [27], we calculate the time derivative of the energy of electron j.
To do so, we apply the chain rule in Eq. ( 2) and obtain where we use ṙj = ∂Ej (t) ∂pj and ṗj = − ∂H ∂rj .In Appendix A, we derive each of the terms in the chain rule in Eq. (20).Furthermore, we group together all the terms in Eq. ( 20) that do not depend on Ėi as follows where and The derivatives and are obtained in Eq. (B1) and Eq.(B9) of Appendix B, respectively.These derivatives are obtained in terms of the relative coordinates q k , since we propagate in the regularised coordinates system.We also find that • ṙi where ṙn = pn mn , ṙi = pi mi and ṙj = pj mj .The derivatives of C j,n and V eff in Eq. ( 24) can be found in Eqs.(B3), (B2) and (B10) in Appendix B. These derivatives and all the terms in Eqs. ( 22) and ( 24) are obtained in terms of the relative coordinates q k , since we propagate in the regularised coordinates system.For each electron we obtain an equation similar to Eq. ( 21).Hence, at any time during time propagation, we solve a system of N -N c equations to obtain the derivative in terms of the energy of each electron, so that it does not depend on the derivatives of the other electron energies.We solve this system of linear equations using Cramers rule [45,46].

E. Hamilton's equations of motion
Substituting Eqs. ( 11) and (18) in Eq. ( 7), we find that the Hamiltonian in regularized coordinates is given by where and p, r, are expressed in terms of ρ and q via Eqs.( 18) and (19).Moreover, we set c k (t) = 0 when k corresponds to the relative distance between an electron and a core.This is the case, since in our model the Coulomb potential between an electron and a core is given by the full Coulomb potential.Using Eq. ( 25), we find that Hamilton's equations of motion are given by where We find ∂V k ′ ,n ∂q k to be given by Thus, we obtain for where the derivatives of V eff and C can be found in Appendix B in Eq. (B11) and Eq.(B4), respectively.Note that when N c = 1, Hamilton's equations in Eq. ( 27) reduce to the corresponding equations for an atom [27].

F. Time derivative of the Hamiltonian
During propagation, besides other criteria used that we do not detail here, we check the accuracy of our propagation also by comparing the propagated Hamiltonian with the Hamiltonian obtained by substituting the propagated variables.The time derivative of the Hamiltonian in Eq. ( 25) is given by The partial derivative of the Hamiltonian with respect to time is given as follows

G. Tunneling during propagation
During time propagation, we allow for each bound electron to tunnel at the classical turning points along the axis of the electric field using the Wentzel-Kramers-Brillouin (WKB) approximation [47].For the transmission probability we use the WKB formula for transmission through a potential barrier [47], with V tun (r, t tun ) the potential a bound electron i can tunnel through given by ϵ i is the energy of the electron at the time of tunneling, t tun , and r a and r b are the classical turning points.We find that accounting for tunneling during time propagation is necessary in order to accurately describe phenomena related to enhanced ionization [48][49][50][51][52] during the fragmentation of strongly driven molecules.

H. Definition of quasifree and bound electron
In the ECBB model, the interaction between a pair of electrons where at least one is quasifree is described with the full Coulomb potential.Effective Coulomb potentials are used to describe the interaction between bound electrons.At the start of time propagation, the tunneling electron is considered quasifree and the other two electrons are bound.We decide on the fly, during time propagation, whether an electron is quasifree or bound using a simple set of criteria described briefly below [27].
A quasifree electron i transitions to bound if, following a minimum approach to the cores, the position of the electron along the field axis is influenced more by the cores than the electric field.We assume that the electron is influenced more by the cores if its position along the electric field has at least two extrema of the same kind in a time interval less than half a period of the laser field.The minimum approach to the cores is identified by a maximum in the Coulomb potential of the quasifree electron with the cores.Also, at the end of the laser pulse, if the quasifree electron has negative compensated energy [53], this electron transitions to bound.In our studies, we use a compensated energy of an electron i that includes the effective potentials as well and is given by A bound electron i transitions to quasifree if either of the following two conditions is satisfied: (i) the compensated energy of the bound electron converges to a positive value or (ii) the magnitude of the total Coulomb potential of the electron i with the cores is smaller than a threshold value and it continuously decreases.The criteria are discussed in detail and illustrated in Ref. [27].
In the initial state of HeH + 2 , all three atoms are placed along the z axis.The two hydrogen atoms are at -3.09 a.u. and -1.02 a.u., respectively, and the helium atom is at 1.04 a.u., with the origin of the coordinate system set to be the center of mass of the molecule.We refer to H farther away from He as H left and the one closest to He as H middle, see Fig. 1.We compute the distance between the two hydrogen atoms and the hydrogen and helium atoms using the quantum chemistry package MOL-PRO [54], employing the Hartree-Fock method with the augcc-pV5Z basis set.The Hartree-Fock method overestimates by a small amount the distance between the hydrogen and the helium atoms [55].However, we employ this method for consistency with the Hartree-Fock wave functions that we use in the potential energy terms involved in computing the exit point of the tunnel-ionizing electron [43].All three nuclei are initiated at rest.

Tunnel-ionizing electron
The electric field is along the axis of the linear molecule, i.e., the z axis with a field strength within the below-the-barrier ionization regime.As a result, one electron (electron 1) tunnel ionizes at time t 0 through the field-lowered Coulomb potential.We employ a quantummechanical calculation to compute this ionization rate as described in Ref. [40].We find t 0 , using importance sampling [56] in the time interval [−2τ, 2τ ] where the electric field is nonzero; τ is the full width at half maximum of the pulse duration in intensity.The importance sampling distribution is given by the ionization rate.We assume , at the time t0 when we initialize our system.The origin of the coordinate system is set to be the center of mass of the molecule.Zcm is the distance between the center of mass of the molecule and the middle of the distance between H left and H middle.
that electron 1 exits along the direction of the laser field; for details on the exit point, see Ref. [43].We compute the first ionization energy of HeH + 2 , with MOLPRO and find it equal to 1.02 a.u.The tunnel ionizing electron exits the field-lowered Coulomb barrier with a zero momentum along the direction of the field.The transverse electron momentum is given by a Gaussian distribution.The latter arises from standard tunneling theory [57][58][59] and represents the Gaussian-shaped filter with an intensity-dependent width.

Microcanonical distribution
In the ECBB model we obtain the initial position and momentum of each bound electron i at time t 0 using a microcanonical distribution with an energy where We take the energy of each electron to be equal to −I p,2 , with I p,2 , being the second ionization potential of the molecule under consideration.We note that the potential energy W i of each electron i in Eqs.(36) and (37) depends not only on the coordinates of the bound electron i but also on the coordinates of all other bound electrons.This dependence is due to the function Hence, the microcanonical distributions of all bound electrons are interrelated, unlike the case for atoms [27].For a triatomic molecule, as is HeH + 2 (the molecule considered here), indices 1-3 are reserved for the nuclei, index 4 for the electron that tunnels in the initial state and indices 5 and 6 for the bound electrons.The microcanonical distribution for the two bound electrons is given by where N is a normalisation constant.To set-up the microcanonical distribution we place the origin of our coordinate system in the middle of the distance between the nuclei A and B, see Fig. 2. As discussed at the end of this subsection, once we obtain the initial conditions for the position of the electrons in this coordinate system, we shift the positions with respect to the center of mass of the triatomic molecule.Moreover, λ, µ are the confocal elliptical coordinates defined using two of the nuclei as the foci of the ellipse with λ ∈ [1, ∞) and µ ∈ [−1, 1] and are given by with r a , r b the relative position between the electron and the two nuclei labelled as A and B, as shown in Fig. 2. R ab is the distance between the two nuclei that are used to define the elliptical coordinates.The coordinate ϕ is the angle between the projection of the position vector r i of the bound electron i on the xy plane and the positive x axis.The z axis goes through the two nuclei that define the elliptical coordinates.Hence, the angle ϕ defines the rotation angle around the z axis, for more details see [9].Transforming from Cartesian to elliptical coordinates, we find that the microcanonical distribution has the following form As in our previous derivation of the microcanonical distribution for one bound electron in the presence of three nuclei [9], we find that the distribution ρ in Eq. ( 41) becomes infinite only when an electron is placed on top of the third nucleus, labelled as C. Hence, as for our derivation in Ref. [9], to eliminate this singularity we introduce an additional transformation with R ac being the distance between nuclei A and C, and similarly for the other distances.We select the value γ = 3 as the lowest one which eliminates the singularity mentioned above [9].The final form of the microcanonical distribution is ρ (λ 5 , t 5 , ϕ 5 , λ 6 , t 6 , ϕ 6 ) = where with where, The parameters x c , z c are given below The new distribution ρ goes to zero when one of the electrons is placed on top of the nucleus of the nucleus C, i.e., when Next, we generate initial conditions for the linear triatomic molecule HeH + 2 assuming the nuclei A,B,C correspond to the Hydrogen atom on the left, to the Hydrogen atom in the middle, and the Helium atom, as shown in Fig. 1.We now identify the range of values of λ i , t i , ϕ i so that P i ≥ 0 for each bound electron i.We find that t i ∈ [t min , t max ], ϕ i ∈ [0, 2π] for each electron, with t min = −(1 + µ c ) 1/γ and t max = (1 − µ c ) 1/γ .That is, P i ≥ 0 is satisfied for the whole range of values of ϕ i and t i .In addition, for P i ≥ 0 to be satisfied we find that λ i cannot be larger than λ max , i.e., λ i ∈ [1, λ max ].The value λ max is the same for both bound electrons.For these range of values, then, we find the maximum value ρmax of the microcanonical distribution ρ (λ 5 , t 5 , ϕ 5 , λ 6 , t 6 , ϕ 6 ) given in Eq. (43).Next, we generate the uniform random numbers λ i ∈ [1, λ max ], t i ∈ [t min , t max ], ϕ i ∈ [0, 2π] for each electron, and χ ∈ [0, ρmax ].If ρ (λ 5 , t 5 , ϕ 5 , λ 6 , t 6 , ϕ 6 ) > χ then the generated values of λ i , t i and ϕ i are accepted as initial conditions, otherwise, they are rejected and the sampling process starts again.
Once we find the λ i , t i and ϕ i , we obtain the position vector r i = (r x,i , r y,i , r z,i ) and the momentum vector p i = (p x,i , p y,i , p z,i ) of each electron i as follows where ϕ p,i ∈ [0, 2π] and ν p,i ∈ [−1, 1] define the momentum p i in spherical coordinates.Following the above described formulation, we obtain the initial conditions of the electron with respect to the origin of the coordinate system.However, for our computations we need to obtain the initial conditions for the position of the electron with respect to the center of mass of the triatomic molecule.To do so, we shift the coordinates by R cm = (X cm , 0, Z cm ) with For HeH 2 , m a , m b are the masses of the hydrogen atoms and m c the mass of the helium atom.We also note that since HeH + 2 is a linear molecule, x c is zero in Eqs. ( 45), ( 46) and (55).Z cm can be seen in Fig. 2. Also, we use the parameters, I p,2 = −1.73a.u., R ab = |r 1 −r 2 | = 2.07 a.u., R bc = |r 2 − r 3 | = 2.06 a.u., R ac = |r 1 − r 3 | = 4.12 a.u., Q 1 = Q 2 = 1 a.u. and Q 3 = 2 a.u.The distances and the second ionization energy of HeH + 2 were obtained using the quantum chemistry package MOLPRO [54], with the Hartree-Fock method employing the aug-cc-pV5Z basis set.

III. RESULTS
Here, we employ a vector potential of the form where k = ω/c is the wave number of the laser field and τ is the full width at half maximum of the pulse duration in intensity.The direction of both the vector potential and the electric field is along the z axis.We take the propagation direction of the laser field to be along the y axis and hence the magnetic field points along the x axis.The intensity of the field is 2 × 10 14 W/cm 2 with a pulse duration in intensity of τ = 40 fs at 800 nm.
Using the ECBB model for molecules, we focus on triple ionization (TI), frustrated triple ionization (FTI), double ionization (DI), and frustrated double ionization (FDI).Out of all events, we find that TI events account roughly for 1.2% and FTI events with n > 2 account for 0.3 %, while DI and FDI with n > 2 events account for 54% and 9.5 %, respectively.Hence, FDI is a major ionization process in strongly driven molecules.In triple ionization, three electrons escape and He 2+ and two H + ions are formed.In frustrated triple ionization, two electrons escape and one electron finally remains bound at a Rydberg state either at He 2+ or at one of the two H + ions.We also find that the formation of He + * + 2H + is roughly 2.5 times more likely than the formation of He 2+ + H + + H * .
In double ionization, two electrons escape, while one remains bound.For the vast majority of DI events, we find that the bound electron has principal quantum number n = 1.Also, it is three times more likely for the final fragments to be He + + 2H + rather than He 2+ + H + + H.That is, in DI, it is three times more likely for the electron to remain bound at He 2+ .In frustrated double ionization, an electron escapes, an electron remains bound at an n = 1 state, and another electron remains bound at a Rydberg state.In FDI, we have several possibilities for the formation of different ions depending on which ions the bound electrons are attached to.We find that it is roughly four times more likely for the deeply bound electron to remain bound at He 2+ versus at H + .
Here, in FTI and FDI we do not include the Rydberg n = 2 states, since an electron from the n = 1 state of H + tunnels to the n = 2 state of He 2+ resulting in a large number of n = 2 states.These states were also not included in our previous work on the strongly driven heteronuclear molecules HeH + [60] and HeH + 2 [40].Finally, we identify the principal quantum number n for each Rydberg electron by first calculating the classical principal quantum number with ϵ i (t f ) being the energy of a bound electron at the end of the time propagation.Then, we assign a quantum number n so that the following criterion is satisfied [61], A. KER distributions In Fig. 3, we plot the kinetic energy release (KER) of the final ion fragments for triple and the most proba-ble route to frustrated triple ionization.As expected, we find that the KER for TI and FTI are very similar.This is consistent with the Rydberg electron in FTI remaining bound in a highly excited state.Hence, the Rydberg electron does not significantly screen the core it remains bound to.Given the similarity of the KER distributions for FTI and TI, for simplicity, we next focus on describing the features of the KER for TI.As in our previous work [40], we find that the left H + ion is the fastest one, followed by He + and the middle H + ion.This is consistent with both Coulomb forces exerted on the left H + ion being along the −z axis.Similarly, both repulsive forces acting on He + are along the +z axis.However, the mass of He is four time larger than the mass of H.This results in He + having a smaller acceleration and hence kinetic energy compared to the left H + ion.Also, the middle H + ion experiences repulsive forces from the left H + ion and from the He + ion in opposite directions.This small net Coulomb force on the middle H + ion results in its kinetic energy being smaller compared to the other two ions.In Fig. 4, we plot the kinetic energy release of the final ion fragments for double and frustrated double ionization.In what follows, we focus on the most probable channels of FDI, namely, the He + + H + + H * [see Fig. 4(b)] which accounts for 65% of FDI and the channel He + * + H + + H [see Figs.4(d) and 4(f)] which account for 17% of FDI.We plot the KER when the bound n = 1 electron is attached to the He 2+ ion (top row), to the left H + ion (middle row), and to the middle H + ion (bottom row).As for TI and FTI, we find that the KER distributions for DI and FDI are very similar.Hence, for simplicity, we next focus on describing the features of the KER distribution for DI.When the deeply bound electron is attached to He 2+ (top row), we find that the left H + is the fastest ion, followed by He + and the middle H + ion.Indeed, the Coulomb repulsive forces on the left H + ion from the other two ions are both along the −z axis, on He + both forces are along the +z axis, and on the middle H + ion the two forces are in opposite directions.When the deeply bound electron is attached to the left H + ion (middle row), the electron screens the charge of the core, resulting in a smaller Coulomb repulsion between the left H fragment and He 2+ .Hence, the kinetic energy of each of the two fragments is smaller compared to their kinetic energy when the deeply bound electron is attached to He + .The reduction in kinetic energy is larger for the left H fragment since both the Coulomb forces from the other two ions are now smaller.This is clearly seen by comparing the grey and light blue lines in Fig. 4(b) with the ones in Fig. 4(a).In contrast, the total Coulomb repulsion on the middle H + ion is increased.Indeed, the repulsion from the left H fragment towards the z axis is decreased while the repulsion from He 2+ towards the −z axis is increased (the n = 1 electron is no longer attached to He 2+ ).Hence, the kinetic energy of the middle H + is increased, compare the purple line in Fig. 4(b) with the one in Fig. 4(a).When the deeply bound electron is attached to the middle H + ion, the repulsive forces on this fragment are smaller resulting in its smaller kinetic energy, compare the purple line in Fig. 4(e) with the one in Fig. 4(c).Also, the kinetic energy of the left H + ion increases since the deeply bound electron now screens the middle H fragment, compare the light blue line in Fig. 4(e) with the one in Fig. 4(c).Finally, the kinetic energy of He 2+ is smaller since the deeply bound electron is attached to the H + ion that is closer to He 2+ , compare the grey lines in Fig. 4(e) and Fig. 4(c) Moreover, we find that the KER distribution of the left H + ion has a pronounced double peak structure for DI and FDI when the deeply bound electron is attached to He 2+ [see light blue lines in Figs.4(a) and 4(b)].To identify the origin of this double peak, it suffices to focus on the DI process.In Fig. 5, we plot the final z component of the velocity of the ions along the electric field, v z , for TI (top row) and DI (bottom row).Also, we plot the contribution to this velocity from the electric field, ∆v E z , and from the forces due to the Coulomb, ∆v C z , and the effective, ∆v V eff z , potentials.Both for TI and DI, we find that the velocity v z of the left H + ion is mostly determined by the Coulomb forces acting on this ion.For DI, the contribution of the Coulomb forces to the velocity of the left H + ion has a clear double peak structure, see Fig. 5(d).This gives rise to the double peak structure in the kinetic energy of the left H + ion in Fig. 4(a).In contrast, the velocities v z of the middle H + and of the He 2+ (He + ) ion for TI (DI) are determined by the forces from both the Coulomb and the effective potentials.The role of the effective potential is more pronounced for the DI process, consistent with one electron remaining bound.Next, for DI, we show how the Coulomb forces exerted on the left H + ion, when the deeply bound electron is attached at He 2+ , result in a low and high energy peak in the KER of the left H + ion, see light blue line in Fig. 4(a).To do so, we plot the angle of escape of each ion with respect to the z axis.The angular distributions that correspond to the low and high energy peaks [see Fig. 6(a)] are shown in Fig. 6(b) and Fig. 6(c), respectively.We find that the two peaks in the kinetic energy distribution of the left H + ion are associated with a different range of angles of escape of the middle H + as well as of the He + ions.For both peaks, the left H + ion escapes along the −z axis, see dark grey lines in Figs.6(b)-(c).Moreover, all ions have roughly the same charge and the middle H + ion is closer compared to He + to the left H + ion.Hence, it is the Coulomb repulsion between the two H + ions that mostly determines the final kinetic energy of the left H + ion.We find that for the lower energy peak the middle H + ion escapes with a very wide range of angles away from the -z axis, compared to a much smaller range for the high energy peak.When the middle H + ion escapes with larger angles with respect to the -z axis and, hence, with respect to the left H + , the Coulomb repulsion between the two H + ions is smaller resulting in a smaller kinetic energy of the left H + ion.In Fig. 7(b), we show that for TI the distributions of the angles of escape of the ions are very similar to the angles corresponding to the high energy peak for DI, compare Fig. 7(b) with Fig. 6(c).As a result, the kinetic energy distribution of the left H + ion for TI is similar to the part of the distribution for DI that corresponds to the high energy peak.Finally, a comparison of the KER for TI, FTI and DI obtained with the ECBB model and with the previous model in Ref. [40] reveals that the KER have larger values for the ECBB model.This is consistent with taking into account the repulsion between the bound electrons using effective potentials in the ECBB model.In our previous more primitive model, the repulsion between bound electrons is turned off.Due to this repulsion via the effective potentials the electrons are less bound to the nuclei they are attached to.Hence, the electrons screen the nuclei less, leading to higher Coulomb repulsion between the nuclei and, therefore, to larger values of the KER.Other differences in the KER of the left and middle H + ions when using the ECBB model versus our previous model in Ref. [40] are due to the effective potentials in the ECBB model significantly influencing the final velocities of the middle H + and He + ions, see Fig. 5.

B. Correlation in electron escape
In Fig. 8, we plot the distribution of the difference of the ionization times between the fastest and second fastest electrons as well as the fastest and slowest electrons in TI and between the fastest and slowest electrons in FTI and DI.We find that the electron that ionizes second has a significant probability to do so with a small time difference from the fastest one, with the time difference being the smallest in TI, followed by DI and then by FTI.The distributions in all three processes extend up to 10 periods (T) of the laser field.In contrast, compared to the fastest electron, the time the last electron ionizes in TI has a distribution that peaks roughly around three periods of the laser field.This suggests that the last to ionize electron escapes mainly due to enhanced ionization and not due to a recollision, i.e. the electronic correlation is weak.This is consistent with HeH + 2 being driven by a long and intense laser pulse in the current study.
A comparison between the distributions of the difference of the ionization times for TI, FTI and DI of HeH + 2 obtained with the ECBB model versus its predecessor model in Ref. [40], reveals that these distributions for the ECBB model peak at smaller times and are less wide.This is consistent with the interaction between bound electrons being accounted for via effective potentials in the ECBB model.As a result, following a return of an electron to the core, energy transfer between bound electrons can take place, leading to possible ionization or excitation.This in turn leads to electrons escaping faster which is consistent with the KER of the fragments having larger values for the ECBB model, as discussed above.

C. n quantum numbers
Next, we investigate the distribution of the principal n quantum number of the two main pathways A and B of FTI and FDI, see Fig. 9.We find that FDI with n > 2 is a major ionization process accounting for roughly 9.5% of all events while FTI is an order of magnitude less probable.Also, we find that pathways A (54 %) and B (46 %) with n > 2 contribute roughly the same to FTI, while for FDI, pathway B contributes significantly more (70%) than pathway A. This can be understood in terms of electronic correlation.FTI, most likely, occurs when the slowest electron finally remains bound in an excited state.As for the slowest electron in TI, see black line in Fig. 8, the slowest electron in FTI can gain energy both from the initially tunneling electron returning to the molecular ion as well as from an enhanced ionization process.Hence, when this electron remains bound in a Rydberg state, it does so either through pathway B, related to energy gain from the returning electron, and through pathway A, related to energy gain from an enhanced ionization process [6].FDI, most likely, occurs when the slowest out of the two electrons that ionize in DI finally remains in an excited state.However, the slowest electron in DI mainly gains energy from the electron returning to the molecular ion, associated with pathway B of FDI.This is supported by the significantly faster ionization time of the second electron in DI and TI (dark and light grey lines in Fig. 8) compared to the ionization times of the slowest electron in TI (black line in Fig. 8).
For FDI, we find that it is significantly more likely for the Rydberg electron to be attached to the H + ion versus the He 2+ ion.Indeed, in Fig. 9(a) and Fig. 9(b), it is clearly seen that the probability for the Rydberg electron to be attached to He 2+ (area under the light grey lines) is much smaller than the probability to be attached to one of the H + ions (area under the dark grey lines).This is consistent with 65% of FDI events having the Rydberg electron attached to one of the H + ions, while the deeply bound electron is attached to the He 2+ , see Fig. 4(b).Only 21% of FDI events have the Rydberg electron attached to He 2+ , while the deeply bound electron is attached to one of the H + ions, see Fig. 4(d) and Fig. 4(f).The significantly higher probability for the Rydberg electron to be attached to one of the H + ions, is consistent with the bound n = 1 electron staying mostly attached to the He 2+ ion both for DI and FDI.In this case, all nuclei have roughly charge of one.In addition, the bound n = 1 electron repels the Rydberg electron from the He + ion, resulting in the Rydberg electron being more likely to stay bound in one of the two H + ions.The Rydberg electron can also remain bound at He + , a less likely process that we do not show in Figs. 4 and 9.
For FTI, in contrast to FDI, it is roughly 2.5 times more likely for the Rydberg electron to remain attached to He 2+ versus the H + ions, compare the area under the light and dark grey lines in Figs.9(c) and 9(d).This is consistent with an electron being significantly more likely to be attracted and remain bound at He 2+ with charge two versus at H + with charge one.
Also, for FDI and FTI, we find that for both pathways the distribution of the n quantum number peaks around n = 18 when the Rydberg electron is attached to He 2+ compared to the significantly smaller n values when the Rydberg electron is attached to H + .This comes as no surprise.Indeed, we assume that the electron that tunnel ionizes last and remains bound in a Rydberg state has roughly the same energy for attachment at He 2+ or at H + .Given that He 2+ has twice the charge of H + , it follows that a Rydberg n state of H + corresponds to roughly a Rydberg 2n state of He 2+ .In addition, we find that the distribution of the n number peaks at higher values for pathway A of FDI versus FTI.This is consistent with an n = 1 electron remaining bound in FDI resulting in a higher screening of the cores in FDI compared to FTI and hence higher energies of the Rydberg electron for FDI.

IV. CONCLUSIONS
We have developed a general three-dimensional semiclassical model for the study of correlated multielectron escape during fragmentation of molecules driven by intense infrared laser pulses.This model fully accounts for the motion of all electrons and nuclei.Moreover, it is developed in the nondipole approximation, fully accounting for the magnetic field of the laser pulse.This model is a generalization of a model we have previously developed for atoms.In this model, referred to as ECBB model for molecules, the interaction of each quasifree or bound electron with the cores and each quasifree electron with a bound electron is treated exactly, fully accounting for the Coulomb singularity.To avoid artificial autoionization, the interaction of a pair of bound electrons is treated through effective Coulomb potentials.We employ the ECBB model in the context of the linear triatomic HeH + 2 molecule.We focus our studies on triple and double as well as frustrated triple and double ionzation.We find that the sum of the final kinetic energies of all ion fragments are larger when described by the ECBB model versus a predecessor of the ECBB model that does not ac-count for the interaction between bound electrons.This suggests that the interaction between bound electrons allows for a more correlated electron-electron escape which occurs faster, leading to a Coulomb explosion of the nuclei at shorter distances.Finally, the ECBB model allows for the study of frustrated double ionization, a major ionization process.This process was not accessible with our previous models, since it has two bound electrons and the interaction between bound electrons was not accounted for.We expect that the ECBB model for strongly driven molecules will pave the way for previously inaccessible studies of multielectron ionization processes during fragmentation of strongly driven molecules.
The function ρ i,n has the following derivatives with respect to r i , r n , q k and ζ i,n ∂ρ i,n ∂r i = 2ζ i,n r n − r i |r n − r i | ρ i,n = 2ζ i,n q k(n,i) q k(n,i) ρ i,n (B5)  ∂V eff (ζ j ′ ,i , (B11)

FIG. 1 .
FIG.1.Schematic illustration of the molecule under consideration, HeH + 2 , at the time t0 when we initialize our system.The origin of the coordinate system is set to be the center of mass of the molecule.Zcm is the distance between the center of mass of the molecule and the middle of the distance between H left and H middle.

FIG. 2 .
FIG.2.The configuration of the triatomic molecule we use to set-up the microcanonical distribution.The origin of the coordinate system is set to be the middle of the distance between the A and B nuclei.

FIG. 3 .
FIG. 3. Distribution of the sum of the final kinetic energies (black solid lines) of the ions produced in (a) triple ionization, and (b) frustrated triple ionization.The grey lines depict the distribution of the final kinetic energy of the He 2+ ion for TI, and He + * for FTI.The purple (light blue) lines depict the distribution of the final kinetic energy of the middle (left) H + ion for TI and FTI.All distributions are normalized to one.

FIG. 4 .
FIG.4.Distribution of the sum of the final kinetic energies (black solid lines) of the ions produced in (a),(c),(e) double ionization, and (b),(d),(f) frustrated double ionization with the n = 1 electron bound at He 2+ (top row), bound at the left H + ion (middle row) and bound at the middle H + ion (bottom row).The grey lines depict the distribution of the final kinetic energy of the He + ion for DI and FDI (top row), and He 2+ for DI and He + * for FDI (middle and bottom row).The purple (light blue) lines depict the distribution of the final kinetic energy of the middle (left) H + ion fragment for DI and H * or H + for FDI (top row), H + (H) ion for DI and FDI (middle row), and H (H + ) ion for DI and FDI (bottom row).All distributions are normalized to one.

5 FIG. 5 .
FIG. 5. Distribution of the z component of the final velocity of the ions, vz, (solid black lines) and of the change of vz in the time interval [t0, t f ] due to the forces from the electric field ∆v E z (dotted green lines), from the effective potential ∆v V eff z (dotted light blue lines), and from the Coulomb potential ∆v C z (dotted dark grey lines).The top row corresponds to TI and the bottom row DI with the n = 1 electron bound at He 2+ .

FIG. 6 . 1 FIG. 7 .
FIG. 6.(a) Distribution of the final kinetic energy of the left H + ion for DI when the n = 1 electron is bound at He 2+ .Angular distributions of the left H + , middle H + and He + ions for DI when the kinetic energy of the left H + is low (b) and high (c).

8 FIG. 8 .
FIG.8.Distribution of the difference of the ionization times between the fastest and second fastest electrons as well as the fastest and slowest electrons in TI and between the fastest and slowest electrons in FTI and DI.

FIG. 9 .
FIG.9.Distribution of the principal n quantum number for pathways A (a) and B (b) of FDI and for pathways A (c) and B (d) of FTI.The distribution of the n quantum number is also plotted separately when the Rydberg electron remains attached to He 2+ for FTI and He + for FDI (light gray lines) and when it remains attached to H + (dark gray lines).