Simultaneous double ionization of C60 molecules in single-cycle pulses

We theoretically investigate simultaneous double ionization of C60 Buckminsterfullerene clusters within the strong field approximation, taking into account two-body effects like Coulomb blocking. Our analysis suggests that for infrared single-cycle pulses, simultaneous double ionization becomes comparable in magnitude to sequential double ionization. Additionally, estimates show that Coulomb blocking weakens with increasing cluster size and field strength.


Introduction
Strong laser field physics underpins various research areas, such as attosecond science (Hentschel et al 2001, Krausz andIvanov 2009), laser machining of micro-and nano-scale structures (Gattass and Mazur 2008, Yanik et al 2004, Farsari and Chichkov 2009, the generation of ultrashort electron pulses (Malka et al 2008, Marceau et al 2013, Krüger et al 2011, Ehberger et al 2015, Zherebtsov et al 2011 for imaging (Sciaini and Miller 2011), and light-driven electronic nanocircuits (Krüger et al 2011, Zherebtsov et al 2011. The strong field process central to all of the above areas is ionization (Keldysh 1965, Delone andKrainov 1991). To a good approximation, ionization in atoms and small molecules can be treated in the single active electron approximation. Electrons are ionized sequentially and do not interact with each other during ionization and subsequent propagation in the continuum. Two notable exceptions have been identified so far. The first multi-electron process is non-sequential tunnel ionization (Fittinghoff et al 1992, Weber et al 2000. The tunneled electron quivers in the continuum and, on recollision with the parent system, kicks out another bound electron. The second process is shake-up/shake-off ionization (Becker andFaisal 2002, Litvinyuk et al 2005), where the tunneling electron shares energy with the remaining electron core and excites or ionizes another electron. This can either occur during a separate Coulomb interaction after the ionization of the first electron or as a simultaneous two-electron process.
If, in the latter case, both electrons are fully ionized in the process, this is also referred to as simultaneous double ionization (SIDI) in strong laser fields. In atoms, theoretical studies have found it to be negligible (Bauer and Rzazewski 1996). However, with the transition of strong field physics from atoms to more complex systems (Krüger et al 2011, Zherebtsov et al 2011, Bhardwaj et al 2003 this question has to be reconsidered. In this work, we theoretically analyze the potential of SIDI in C 60 , representative of complex systems. Our analysis is based on the strong field approximation (SFA), in which the ionic potential of the parent ion is neglected. We find that for infrared single-cycle pulses, SIDI becomes comparable in magnitude to sequential double ionization (SEDI). Finally, we find qualitatively that the relative difference between SIDI and SEDI shrinks even more with increasing field strength and cluster size. This is reasonable as two electrons in the ground states are further apart and affect each other less during tunneling.

Theoretical approach
We use atomic units throughout the paper unless stated otherwise.

Correlated simultaneous double ionization
For SIDI, we consider two active electrons moving in the mean field created by the non-active electrons and the atomic cores. Using non-relativistic dipole approximation and length gauge, the Hamiltonian of the system reads where F(t) is the time-dependent electric field in dipole approximation, t i is the kinetic energy operator of the ith electron, v i is the effective potential energy operator accounting for the interaction with the atomic cores and the non-active electrons, and w(1, 2) = 1/|r 1 − r 2 | is the interaction energy operator between the active electrons. An important approximation in equation (1) is the neglect of the cluster polarization field which would ensure that the total field strength at the cluster surface is zero. For a given cluster, this approximation is justified as long as the tunneling barrier is not suppressed below the binding energy of the ground state.

Initial state-rigid-rotor model
The ansatz for the 2-electron ground state present before the onset of the electric field is obtained using the rigid-rotor approximation (Savina et al 1993) that describes the electrons of C 60 as a free 2D electron gas on a sphere of radius r 0 . Using the symmetry adaption to the icosahedral symmetry of C 60 (Qiu and Ceulemans 2002), the single-electron states are superpositions of the form where are the rigid-rotor eigenstates of a sphere with radius r 0 , are the spherical harmonics with the normalization factor N lm , and P l|m| (·) the associated Legendre polynomials. u is the quantum number of the symmetry-adapted state, which is a superposition of rigid-rotor states of different m with real-valued coefficients p um . These coefficients are given explicitly e.g. in (Qiu and Ceulemans 2002). If not stated otherwise, r 0 ≈ 3.5 Å (Goel et al 2004) is used as the radius of the C 60 molecule throughout this work. The initial state present at a time t i before the pulse is approximated by a single 2-electron configuration of rigid-rotor states, i.e.
where |s i , m s i is a coupled 2-electron spin state, and is a normalization factor. As the HOMO level electrons' angular momentum in the rigid rotor model is l = 5 (Savina et al 1993), we are using l 1 = l 2 = 5 throughout this work. Thus, i = {u 1 , u 2 , s i , m s i } denotes the set of quantum numbers fully describing the initial state of the 2-electron system. We assume that where the ground state energy E d ≈ −18.98 eV is the energy of the initial 2-electron state, i.e. the opposite of the double ionization potential (Pogulay et al 2004) between the ground state and the 2-electron continuum at E = 0.

Final state-correlated Coulomb wave
In order to fully include correlation in the final state, an eigenfunction of the Hamiltonian is used as the final state present after the pulse at time t f . This Hamiltonian describes 2 free electrons interacting with each other due to their Coulomb repulsion. Its eigenstates read (Landau and Lifschitz 1979, pp 101-129) in position-space representation using center-of-mass (COM) and relative coordinates R and r = (r, θ, φ), respectively. Here, is the real-valued Coulomb wave function with the normalization factor where Γ(·) is the Gamma function and F (·, ·, ·) is the confluent hypergeometric function. Due to the fermionic antisymmetry of the final state under particle exchange, l + s f is always even. The corresponding eigenenergies read where M and μ are the total and reduced masses of the 2-electron system, respectively. Note that the parameter k includes both the kinetic and the potential energy of the relative motion of both electrons. As the ionized electrons will move away from each other for t → ∞, and since after ionization, the relative and COM coordinates of the 2-electron system do not interact anymore, the potential energy of the relative motion will eventually be entirely converted to kinetic energy. Therefore, k can be interpreted as the asymptotic momentum of the relative motion. We denote the set of discrete and continuous quantum numbers fully describing the 2-electron system's final state with f = {K, k, l, m} in the following.

SFA
In SFA (Amini et al 2019, Keldysh 1965, Faisal 1973, Reiss 1980, it is effectively assumed that the exact shape of the molecular binding potential is negligible. Hence, only the value of the ionization potential is included in the calculation. Using SFA in length gauge, the ionization amplitude for SIDI M SIDI if in a time-dependent electric field F(t) in non-relativistic dipole approximation reads where is the dipole moment, is the classical action, and is the vector potential of the electromagnetic field. Note that the factor 2 in front of the vector potential A in (13) and (15) as well as in front of the electric field F in (13) is different from single ionization and is due to the fact that the 2-electron system is doubly charged.
We consider a linearly polarized single cycle pulse. More precisely, we assume a field of the form The model for SIDI that we introduce below is valid in the quasi-static limit. Thus, it is assumed that in a multi-cycle pulse, the ionization of the first half cycle calculated using our model is repeated into the opposite direction during the next half-cycle, and then the process repeats for the following cycles.

Integration using the saddle-point approximation
The integral over t in (13) can be simplified further using the method of steepest descent (saddle point approximation). Here, we assume quasi-static dynamics, i.e. a small Keldysh parameter (Keldysh 1965) For atomic systems, it is typically assumed further that F 0 Z |E d | for all values of Z that are covered by the initial state, from which would follow that the term oscillates much slower than exp(−iS f (t )) in the integration variable t . For the type of systems this work is concerned with, i.e. large molecular systems like C 60 in strong fields, this assumption would however not be correct. Instead, we break down the transition dipole moment ψ lm (K, k)|Z|0 i into its longitudinal Fourier componentsd if (Z) according to and then apply the saddle point approximation to each Fourier component separately. This yields is the curvature of the exponent at the saddle point, and is a phase ensuring that the integration path crosses the saddle point in the direction of steepest descent. Note that due to the fact that the dipole-approximated field only couples to the COM coordinates, this expression resembles the tunneling amplitude of a single particle of mass M and charge 2e into a free continuum state with momentum K. The important difference however is that the correlation in the final state is still hidden in the Fourier componentsd if (Z) of the transition dipole moment, and in order to obtain M SIDI if , a six-dimensional integral over all 2 × 3 space dimensions has to be calculated. To do this efficiently, we use three different approximations.
First, since the exponential term in (21) declines exponentially for transversal COM momenta K ⊥ = 0, the final wave packet's (see (9)) transversal COM momentum is very narrowly distributed around K ⊥ = 0, and we can write exp(iKR) ≈ exp(iK Z Z).
Second, the exponential expression in (21) is pronounced for values of the COM position Z that are close to the minimum value with nonzero probability in the initial state. As, within the rigid-rotor model, both electrons are confined to a sphere of radius r 0 , the minimum COM position with nonzero probability in the initial state is Z = −r 0 . Since the expressions Γ f (Z) and Ξ f (Z) only depend weakly on Z, we can Third, we apply an approximation to the radial component of the final state, replacing the radial coordinate r with an effective value The value r eff f describes the most probable separation of both electrons immediately after tunneling and is a function of the final state's quantum numbers. A detailed justification and discussion of this approximation is given in the appendix A. r eff f is shown in figure 1 for the laser parameters F 0 = 0.04 a.u. and λ = 8000 nm. The effective separation r eff f of both electrons decreases with increasing k and l, corresponding to a higher energy of the relative motion, and it is approximately independent of the final state's COM momentum K. Thus, it is assumed that r eff f ≈ r eff kl , where K perp = 0 and K Z = −2F 0 /ω fixed. Using these three approximations in the evaluation of integral (21), we obtain dx 2 (x 1 + x 2 ) P l ,|m 1 | (x 1 )P l−l ,|m 2 | (x 2 )P l 1 ,|m 1 | (x 1 )P l 2 ,|m 2 | (x 2 ) where a numerical factor accounting for the shift of the spherical harmonics from the COM frame to absolute coordinates. The remaining two-dimensional integral over x 1 and x 2 can be evaluated numerically.

Uncorrelated sequential double ionization
In the numerical experiments shown in the next section, we compare photoelectron momentum spectra after correlated SIDI to those after uncorrelated SEDI. For SEDI, it is assumed that the two electrons ionize sequentially and do not influence each other. We emphasize however that for infrared single-cycle pulses as the one used in this work, the analysis in (Bart et al 2020) indicates that Coulomb blocking (CB) effects can play an important role for SEDI of C 60 as well, suppressing SEDI by two to three orders of magnitude. In that sense, the SEDI rates obtained here must be seen as an upper estimate. Starting from the same initial state |0 i , the ionization of each electron is described independently of the other electron, neglecting any Coulomb interaction of the electrons. The single ionization steps are calculated in SFA, in an approach analogous to the one presented for SIDI.
The single ionization amplitude for the jth electron from an initial rigid-rotor state |l j , u j into a plane wave state |k j with wave vector k j reads where denotes a real-valued phase, t b,j = ω −1 arccos(−1 − ωk z,j /F 0 ) is the jth electron's birth time and |E s,j | is the single ionization potential of the jth electron.
is the curvature at the saddle point, and ξ k j (z) = exp(−i/2(arg(γ k j (z)) + π)) is a phase ensuring that the integration path crosses the saddle point in the direction of steepest descent. Note that the amplitude itself has an index j, since the single ionization potential |E s,j | differs between the first (j = 1) and the second (j = 2) ionization process. The total ionization amplitude for both steps reads where f = {k 1 , k 2 } are the momenta at the detector of the first and second electrons to ionize, respectively. As both electrons are accelerated in negative z-direction by the field after tunneling, the second electron's momentum is always larger ('less negative') at the detector. The Heaviside function Θ(·) accounts for this fact. Please keep in mind that speaking of a 'second electron' here is slightly inaccurate, since of course both electrons are interchangeable. The 'second electron' is a term used to describe those parts of the wave function that ionize during the second ionization step with ionization potential |E s,2 |. In the following, the amplitude M SEDI if will be expressed using COM K = k 1 + k 2 and relative k = 1 2 (k 2 − k 1 ) momenta.  Figure 2 shows a comparison of the differential ionization probability

Differential ionization probabilities
for SEDI (figures 2(a) and (c)) to the differential ionization probability for SIDI (figures 2(b) and (d)). We would like to highlight two important distinctions between the electron momentum distributions of SEDI and SIDI that could potentially be used to identify SIDI in an experiment. First, the most likely birth time of the two-electron wave packet in SIDI is at the peak of the laser pulse. Therefore, the distributions resulting from SIDI are symmetric around K z = 2F 0 /ω, which corresponds to a birth time at this peak. In SEDI, the most likely outcome is that the first electron is born before and the second after the peak, meaning that the photoelectron distribution exhibits no such symmetry. Second, the relative momentum k of both electrons depends on the wavelength λ for SEDI, since longer time differences between first and second tunneling are possible for long wavelengths. For SIDI, the relative momentum k is to a good approximation independent of the wavelength. Instead, the distribution in k is determined by an interplay of the size of the initial state and the interaction energy in the final state. We discuss this CB effect in more detail in section 3.2. For this reason, the momentum distributions for SIDI are unchanged along the k-axis (notice however that the K z -axis is scaled with λ −1 ) between λ = 1000 nm  (33) together with their product, f result (k). The laser parameters are F 0 = 0.04 a.u. and λ = 8000 nm. and λ = 8000 nm. For SEDI, this is not true, and much higher vales of k occur with significant probability at longer wavelengths.

The CB effect during SIDI
Small values of k in the final state correspond to a final state with little energy stored in the relative motion of both electrons. This corresponds to a final state in which both electrons are far apart from each other after SIDI, since such a state has low potential Coulomb energy. Conversely, large values of k correspond to both electrons being very close to each other, with a high Coulomb interaction energy. The ground state is confined to a sphere of radius r 0 , which puts an upper limit to the separation of both electrons before and right after tunneling. This translates into a lower limit of the values of k in the final state. On the other hand, large values of k in the final state corresponding to very small separations are also prohibited, as this increases the final state's energy and thus heightens the tunneling barrier. The interplay of both these effects-the initial state's confinement pushing towards larger values of k and the energy term pushing towards smaller values of k-results in a much sharper peak of the distribution in k for SIDI compared to SEDI at λ = 8000 nm.
Mathematically, the initial state's confinement is reflected in the analytical expression for M SIDI if by the remainder of the Coulomb wave function with the effective radius, R kl (r eff kl ). The tunneling barrier prohibiting large values of k manifests in the tunneling exponent, of which is the predominant contribution. Figure 3 shows R k, l=0 (r eff k,l=0 ) and f E (k) together with their product, f result (k). Both factors grow (decline) exponentially for small (large) values of k respectively, thus the resulting product's peak in k is very sharp. Also note that the peak in k coincides approximately with the one in the full distribution in figure 2(b). Since the confinement of the initial state 'pushes' the peak in k towards higher values, the interaction energy of the final state is higher compared to the case in which this confinement is not present. This additional energy suppresses the ionization probability, and it is therefore a manifestation of the CB effect. Furthermore, the strength of this suppression, i.e. the strength of the CB effect during SIDI, is dependent upon the extent to which the initial state is confined, i.e. upon the geometry and the size of the initial ground state.

Varied laser parameters
We compare SIDI to SEDI at varied wavelengths λ and peak laser field strengths F 0 . Figure 4(a) shows the ratio of the total ionization probabilities P SIDI /P SEDI , and figure 4(b) shows the ratio of the peak values of the differential ionization probabilities d 4 P SIDI d 3 Kdk K ⊥ =0 and d 4 P SEDI . Our model is applicable to laser parameters in the quasi-static regime, in which γ < 1. For the field strength F 0 = 0.04 a.u., the Keldysh parameter assumes the value γ = 1 at λ ≈ 925 nm, at F 0 = 0.02 a.u., this border is at λ ≈ 1850 nm. . The second x-axis indicates the Keldysh parameter γ for F 0 = 0.04 a.u. Results for γ 1, to which our model is not rigorously applicable, are indicated as gray crosses. All results are calculated for the initial configuration i = {u 1 = ζ, u 2 = ζ, s i = 0, m s i = 0}. In the case of SIDI, all symmetry-allowed l-contributions up to l max = 12 have been considered as the final state.
The SEDI channel is larger in magnitude than the SIDI channel; the ratio of the peak values shown in figure 2 is ca. 0.016 at λ = 8000 nm and F 0 = 0.04 a.u. At the same laser parameters, the difference in the integrated ionization probabilities P SIDI and P SEDI is around three to four orders of magnitude. While the peak ratio is approximately independent of the wavelength, shorter wavelengths narrow the time window between the two sequential tunneling events in SEDI. Therefore, the relative importance of SIDI increases for shorter wavelengths; at λ = 1000 nm and F 0 = 0.04 a.u., the difference between P SEDI and P SIDI is reduced to two to three orders of magnitude. Importantly, as mentioned before, the CB effect was neglected in SEDI here, and recent work suggests that the CB effect plays an important role for SEDI in short pulses as well. Specifically, in (Bart et al 2020) it was shown that for sub-nm clusters driven by single-cycle infrared lasers, this suppression can be two to three orders of magnitude. In combination with that, our analysis reveals that for sub-nm clusters driven by single-cycle infrared lasers, one can expect the integrated ionization probability of SIDI to become comparable to SEDI.

Varied cluster size
The geometry dependence of SIDI mentioned in section 3.2 can be seen clearly by repeating the calculation presented in section 2 for artificially chosen different molecule radii r var instead of the actual molecule radius r 0 , while freezing all other system parameters. We compare SIDI to SEDI at varied radii r var and peak laser field strengths F 0 . Figure 5(a) shows the ratio of the total ionization probabilities P SIDI /P SEDI , and figure 5(b) shows the ratio of the peak values of the differential ionization probabilities d 4 P SIDI The magnitude of SIDI depends very sensitively on the molecule's size. The difference between SIDI and SEDI shrinks with growing cluster size, since the two electrons in the ground state are further apart and affect each other less during tunneling. Again, recall that the SEDI values reported in figure 5 do not include CB effects, which would further suppress SEDI compared to SIDI. We have argued in section 3.1 that the integrated ionization probability of SIDI is expected to be comparable to the integrated ionization probability of SEDI for C 60 exposed to single-cycle pulses of wavelength 8000 nm and peak field strength F 0 = 0.04 a.u. already. Figure 5 suggests that one could expect such a regime with even higher confidence at higher field strengths or in larger clusters.
For larger clusters, the gap between the first and second ionization energy narrows. The tunneling exponents of SIDI and SEDI are expected to become comparable, since they fulfill . All results are calculated for the wavelength λ = 8000 nm and the initial configuration In the case of SIDI, all symmetry-allowed l-contributions up to l max = 12 have been considered as the final state.
where the equality sign holds if E s,1 = E s,2 . At the same time, in the limit of large clusters, the CB suppression in the SIDI channel becomes weaker due to the larger distance of both electrons after tunneling. As discussed before, SEDI driven by single-cycle laser pulses is weakened by Coulomb blocking (Bart et al 2020). Therefore, the size-dependent weakening of the CB effect applies to SEDI as well, but plays less of a role here: During SEDI, both electrons tunnel at different times, and thus have an additional mechanism to avoid each other, a mechanism that is independent of the size of the molecule. Beyond this qualitative analysis, our model is however not suitable to analyze this regime quantitatively, since scaling it to even larger cluster sizes and field strengths would result in an unphysical suppression of the tunneling barrier. As was discussed in the beginning, this suppression is due to the neglect of the cluster polarization field which makes the total electric field zero at the cluster surface. This effect has not been considered in our analysis.

Conclusion
We have presented an analytical model for SIDI in C60 clusters. Our analysis suggests that SIDI becomes comparable to SEDI in C60. SIDI further gains in importance when the size of the cluster is increased, and therefore is expected to become observable for cluster sizes of a few angstrom and larger. To investigate clusters much larger than C60, polarization effects need to be included into our model. This will be subject to further research. Finally, SIDI specific signatures have been discussed which could be used for its measurement. Figure A2. Comparison of L f (x min (r), x min (r), r) to L f (x min (r), x min (r), r eff kl ), shown for several values of final-state relative asymptotic momentum k and final-state angular momentum l. The final-state COM momentum is fixed at K = (0, 0, −2F 0 /ω). The laser parameters are F 0 = 0.04 a.u. and λ = 8000 nm. denotes the dominant contribution in r to the final three-dimensional integral, and in our approximation, it is only dependent on the parameters f of the final state. Figure 1 shows the value of r eff kl as defined in (A.10) for varying final-state parameters l, k, K ⊥ and K Z for the peak field strength F 0 = 0.04 a.u. The dependency on the final-state parameters K ⊥ and K Z is relatively small and will be neglected, thus r eff f ≈ r eff kl , where K perp = 0 and K Z = −2F 0 /ω fixed. dx 2 P l |m | (x 1 )P l−l ,|m−m | (x 2 )P l 1 ,|m 1 | (x 1 )P l 2 ,|m 2 | (x 2 ) The remaining integrals over x 1 and x 2 have to be calculated numerically. Note that the term δ m 2 ,m−m only occurs because of our approximation, which means that the resulting selection rule is approximate as well.
To verify the applicability of our approximation r 0 r l R kl (r) ≈ r 0 r eff kl l R kl (r eff kl ), we compare the leading term L f (x min (r), x min (r), r) to its counterpart after the approximation is applied, i.e. L f (x min (r), x min (r), r eff kl ). Both terms are shown in figure A2. For values of k and l contributing significantly to the final state (note the absolute scale), the approximation works fairly well and allows us to greatly simplify the numerics.