Interatomic coulombic decay rate in endohedral complexes

Interatomic coulombic decay (ICD) in van der Waals endohedral complexes was predicted to be anomalously fast. However, the available theoretical calculations of the ICD rates in endohedral complexes only consider the equilibrium geometry, in which the encapsulated atom is located at the centre of the fullerene cage. Here we show analytically that the dominant contribution of the dipole plasmon resonance to ICD does not deviate from its equilibrium geometry value, while contributions of higher multipole plasmons to the ICD can be neglected for most atomic displacements possible for an endohedral complex at room temperature. This is in contrast to the behaviour predicted for ionic endohedral compounds. Our results show that the conclusion of the earlier works on the ultrafast character of the ICD in endohedral complexes holds generally for a wide range of geometries possible under a thermal distribution, rather than only for the equilibrium geometry.


Introduction
Intra-atomic non-radiative decay processes [1], such as Auger effect, are most commonly observed in core-ionised atoms and molecules and typically occur on the time scale of a few femtoseconds. On the other hand, an inner-valence vacancy in an isolated ion, whose energy lies below the double ionisation potential of the ion, may only decay via the slow (nanosecond scale) photon emission process. The situation changes dramatically for the inner-valence vacancy in an ion placed in the environment of a weakly bound cluster, in which case an alter- 4 Authors to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. native fast decay channel called interatomic coulombic decay (ICD) [2] becomes open. In the course of the ICD, recombination of an inner-valence vacancy in a sub-unit belonging to a cluster of atoms/molecules causes an electron in the outervalence of another sub-unit to ionise, resulting in two separated charges. In weakly bonded systems, such as van der Waals or hydrogen-bonded clusters, charge separation between cluster sub-units significantly lowers the double ionisation threshold of the cluster in comparison to the isolated atoms or molecules, making ICD energetically possible.
The lifetime of the inner-valence vacancy and the ICD rate depends on the internuclear distance and the number of atoms interacting with the vacancy. For example, the lifetime of the 2s hole on Ne is 85 fs for Ne 2 but only 3 fs for a central Ne atom in Ne 13 [24]. For greater number of neighbouring atoms and hence also decay channels, the ICD rate is higher and the lifetime of the inner-valence vacancy is shorter. This motivates the study of ICD in complex structures such as endohedral fullerenes where the encapsulated atom or molecule interacts with about 60 nearest neighbours, thus maximising the ICD rate. Indeed, initial theoretical study of ICD in endohedral complexes [25] predicted extremely fast ICD, approximately 2 fs for (2s 1 )Ne + @C 60 , which is faster than, for example, intra-atomic Ne 1s Auger decay. It was noted that the energy, released from the recombination of the inner-valence vacancy on the Ne atom, excites the collective mode of the electrons on the fullerene cage known as 'giant plasmon resonance' [26,27]. It is expected that the predicted high ICD decay rate in (2s 1 )Ne + @C 60 has its origin in the near-resonant coupling to the giant plasmon [16,25,28,29]. Theoretical study of ICD in endohedral complexes is further stimulated by the very recent observation of ICD in ionic endohedral compounds [30] (where the endohedral atom forms an ionic bond with the fullerene), which shows that in line with the theoretical predictions, ICD competes successfully with the intra-atomic decay.
The earlier calculations of the ICD rates in van der Waals endohedral complexes were performed at the equilibrium geometry of the cluster, which for van der Waals complexes is at the centre of the fullerene cage [25,28]. Nevertheless, in van der Waals endohedral complexes the confined atom, if it is small enough, moves in a relatively flat potential and can deviate significantly from the equilibrium position. For example at the room temperature, Ne in Ne@C 60 populates a sphere of the radius of about 1 a.u. [31]. Moreover, in ionic endrohedral compounds (where the equilibrium geometry is off-centre) the interatomic decay width has been predicted to vary strongly with the displacement from the central geometry [32]. The question arises, how much is the ICD rate affected by the position of the endohedral atom within the cage? One aspect to consider here is the symmetry breaking. Indeed, as atom moves away from the centre of the fullerene, it could excite not just the dipole plasmon, but also the higher angular momentum modes of collective electronic motion [33]. In order to answer this question, we develop an analytical approach to the ICD rates in van der Waals endohedral fullerene complexes at non-symmetric 'off-centre' geometry.

Theory
We begin by considering the ICD width in the Wigner-Weisskopf approximation without exchange [34] for a single decay channel of the system A + @C 60 , where A is the endohedral atom bearing an inner-valence (iv) vacancy. Here the direct matrix element V direct is where ov 1 (centred on A) and ov 2 (centred on C 60 ) are the two outer-valence orbitals characterising the final state of the decay, while κ is the continuum orbital of the ICD electron. The direct matrix element (2) embodies the physical mechanism of the ICD process: an outer-valence electron (1) recombines into the inner-valence orbital of A, while an outervalence electron (2) is ejected into the continuum as a result of the Coulomb interaction. The corresponding exchange matrix element for this decay channel can be written as Whereas the direct contribution (2) takes place as a result of the long-range effect of the Coulomb force, the exchange contribution (3) requires overlap between the valence orbitals of the two species involved in the process. Here we assume that A is weakly interacting with the fullerene cage and as a result is free to move inside it, as for example in the case of Ne or N [31,35]. This means that there is hardly any orbital overlap between the endohedral atom and the cage for a range of displacements from the central geometry and as a result, the exchange contribution (3) at these geometries can be neglected. In what follows we therefore concentrate on the study of the direct matrix element (2) as a function of the displacement of A from the centre of the fullerene. The direct matrix element (2) can be represented as an electrostatic repulsion of two charge distributions (transition densities) possessing zero total charge each (due to molecular orbital orthogonality). If these distributions were independent of the endohedral atom displacement, the resulting repulsion would have been displacement-invariant as well as it is wellknown in electrostatics. However, while the bound orbitals can be assumed to depend very weakly on the atomic displacement, the continuum orbital of the ICD electron coupled to the displaced bound orbitals by the electron repulsion can change strongly as the endohedral atom departs from the central position (e.g. due to symmetry breaking) and the electrostatic argument strictly does not apply. Thus one has to develop a theory of the displacement dependence of the direct matrix element V direct .
Following earlier analysis of diatomics [34,36] and symmetric endohedrals [28], we would like to express the direct The fullerene cage thickness ΔR is negligible comparing to the average radius (r C 60 ≈ R), R is 6.639 a.u. [38,39]. matrix element V direct in terms of physical quantities that belong to A and C 60 exclusively. To this end, we use multipole expansion of the coulombic repulsion [37]: where electron (1) described by the radius-vector r is assumed to be localised on A, while electron (2) described by the radiusvector r C 60 is localised on the fullerene (see figure 1) in accordance with the zero-overlap assumption. Substituting the multipole expansion (4) into the direct matrix element (2), the ICD matrix element can be written in the form: One can see that the direct ICD matrix element (5) is factorised into two separate contributions characterising the ICD in the following way: the first endohedral matrix element corresponds to the recombination of the inner-valence vacancy on A, while the second (fullerene) one represents the ionisation of C 60 as a sum over multipole transitions. In equation (5) the atomic orbitals on A are expressed as functions of the electron (1) coordinates r with the origin at the fullerene centre. However, analytical expressions for atomic orbitals are given in the coordinate system residing on the atom, in which electron (1) is described by r A , see figure 1. Thus one needs to be able to express the orbitals at the shifted coordinates r A = r − D using the spherical harmonic basis residing at the centre of the fullerene. While we are not aware of a general solution of this problem, it turns out to be possible for Slater-type orbitals (STOs) [40] Ψ ζ n nLM (r A ) = r n−1 A e −ζ n r Y M L (θ, φ), where ζ n = Z eff n is the orbital exponent reflecting the effective charge experienced by the electron [41]. According to Silverstone [40], the spherical-harmonics expansion of the displaced STO can be written as where 2λ + 1 4π (8) are referred to as the Condon-Shortley (or Gaunt) coefficients [37] and functions v ζ n l,λ,L (r, D, n) = (9) carry all the radial dependence of the wavefunctions. Here j λ (kD) and j l (kr) are the modified spherical Bessel functions [40] of orders λ and l respectively, k is the wavenumber, where n is the principal quantum number corresponding to equation (6), and Γ(n − L + 1) is the gamma function. Transition density on the atom A is a product of the ov 1 and iv orbitals each of which can be presented as a linear combination of STOs (6) of the given orbital momentum and its projection. Therefore, the transition density is by itself a sum of STOs: Making an approximation that r C 60 ≈ R in the denominator of equation (5), where R is the radius of C 60 (see figure 1), as the valence carbon electrons are assumed to be confined within the thin layer ΔR of the fullerene cage, and using the representation (11) for the transition density on the endohedral atom, we can rewrite the direct ICD matrix element (5) in the form where the coefficients P lm (D) defined as (13) govern the relative strength of the multipole contributions to the ICD width in A + @C 60 or in other words coefficients P lm (D) play the role of complex amplitudes, the squared absolute value of which gives the probability to excite each plasmon mode corresponding to the angular momentum l and its projection m. Moreover, the multipole weight coefficients (13) alone carry all the dependence of the process on the displacement D of the endohedral atom.

Results
We consider the displacement dependence of the multipole weight coefficients (13) in the typical case, where the innervalence atomic orbital iv is of 2s type, while the outer-valence atomic orbital ov 1 is of 2p type, e.g. 2p z (see figure 1), in the case of Ne atom: here ζ = 2.8792 a.u. [42]. Then the product of STOs (11) entering the weight coefficients (13) for orbitals (14) can be written as iv(r − D)ov 1 D)).
(15) Using the spherical-harmonics expansion of the displaced STOs (see equation (7)) and substituting the result of (15) into the expression for the weight coefficients (13), for the case of Ne (L = 1, M = 0) we obtain: There are two distinct cases for the orbital orientation dictated by the symmetry of the system: (i) Ne's displacement parallel to the z-axis and (ii) perpendicular to the z-axis. Figure 2 presents the behaviour of non-zero relative strength (16) of different modes (for angular momentum values up to l = 4) of the plasmon oscillations excitation as a function of the displacement in the z-direction. One can see that in case of the displacement along the z-direction the oscillations with only m = 0 appear. The main contribution is the dipole-dipole (l = 1) one. This is the only contribution allowed by symmetry at zero displacement and according to the electrostatic (Gauss theorem) argument above, it should not vary with displacement. Numerical evaluation of equation (16)   indeed confirms negligible variation of the dipole-dipole coefficient with D. The contributions drop significantly for higher l values. The monopole l = 0 contribution is not presented as it is cancelled by the zero monopole term on C 60 in equation (13) due to molecular orbital orthogonality.
These features are similar for the results obtained for the situation when the Ne orbital is displaced along the x-axis (also valid for y-axis due to the symmetry of the system) which are shown in figure 3. In contrast to the case of displacements along z-direction, for perpendicular one the other modes of the plasmon oscillations with non-zero m are excited. As in the case of the z-displacement, variation of the dominant dipole-dipole contribution is negligible in compliance with the Gauss theorem argument based on the electrostatic analogy for the direct ICD matrix element (2). This leads to the conclusion that the giant dipole resonance is the only mode of the plasmon oscillation that gets excited following ICD in (2s 1 )Ne + @C 60 when the Ne stays within a spherical region of radius up to 2 a.u. around the centre of the C 60 cage.

Conclusions and outlook
We have shown that the rate of ICD in van der Waals endohedral fullerene complexes is stable with respect to the deviation of the endohedral atom from the fullerene centre. Our analytical theory describes the direct contribution to the ICD width as a weighted product of multipole excitations on the fullerene (plasmons), where the weights define the dependence of the width on the complex geometry. We have shown that for the typical parameters of the inner-valence-ionised complex, e.g. for (2s −1 Ne + C 60 ), the dipole-dipole contribution to the ICD is by far the leading one over the full range of the geometries accessible at room temperature.
The presented multipole expansion analysis of the ICD width in endohedrals is valid under the assumption of zero overlap between the orbitals of the encapsulated atom and the fullerene cage. In cases where this assumption breaks down, e.g. in Ar@C 60 [43], we would expect the ICD to be strongly enhanced relative to the predictions of the present theory, similarly to the behaviour of the ICD rate in the rare gas-alkaline earth diatomics at short distances [34] as is indeed confirmed by the recent calculations [44].
In view of the recent experimental breakthrough on ICD in ionic endohedral compound [30], experimental realisation of the ICD in van der Waals endohedral complexes appears to be within reach. We hope that our work will motivate more active experimental research on ICD in endohedrals, using synchrotron and/or x-ray free-electron laser sources to create the inner-valence endohedral vacancies. ICD could then be detected through the photoelectron-ICD electron [6] or ICD electron-parent molecular ion [30] coincidence spectroscopy. Since dissociative pathways are energetically accessible to the ICD final states [45], a comparative study of the fragmentation of the endohedral complex and the pristine fullerene can serve as an additional indication of the interatomic decay process.