Electronic friction coefficients from the atom-in-jellium model for $Z=1-92$

The break-down of the Born-Oppenheimer approximation is an important topic in chemical dynamics on metal surfaces. In this context, the most frequently used"work-horse"is electronic friction theory, commonly relying on friction coefficients obtained from density functional theory (DFT) calculations from the early 80s based on the atom-in-jellium model. However, results are only available for a limited set of jellium densities and elements ($Z=1-18$). In this work, these calculations are revisited by investigating the corresponding friction coefficients for the entire periodic table ($Z=1-92$). Furthermore, friction coefficients obtained by including the electron density gradient on the Generalized Gradient Approximation (GGA) level are presented. Finally, we show that spin polarization and relativistic effects can have sizeable effects on these friction coefficients for some elements.

Since the BO-approximation is a very fundamental approximation in theoretical chemistry, going beyond imposes a severe conceptual challenge. Alternatively, solving the fully coupled electron-nuclear time-dependent Schrödinger equation to completely avoid the BO approximation altogether will remain computationally intractable for the foreseeable future, even for systems with only very few degrees of freedom. For going beyond the BO approximation, a commonly used approach is combining ab initio molecular dynamics with electronic friction theory 16,20 . Using the local density friction approximation (LDFA) is a way to include the dissipative effect of electron-hole pair excitations in molecular dynamics 14 that is computationally much more convenient than other approaches 21,22,[38][39][40][41] . Within the LDFA including the independent atom approximation, the so-called electronic friction coefficient η is required. η only depends on the nuclear charge of the moving atom and the electron density of the metal surface at its point-like nucleus or different atoms-in-molecule decompositions of the latter 17,42,43 . The friction coefficient is obtained by using the atom-injellium model, where the atom is embedded in an infinitely extended homogeneous electron gas of that density. The energy loss in the jellium model is caused by the momentum loss of the nucleus due to the scattering of the electrons from the gas. The electronic friction coefficient is obtained from the electronic structure of the atom in jellium, which is obtained from density functional theory (DFT) 44,45 .
The local (spin) density approximation (L(S)DA) for the exchange-correlation functional in DFT is by construction exact for the jellium background, which is why the electronic structure of atoms in jellium has traditionally also been obtained at this level of theory 46-48 . However, the electronic structure of the atom in jellium is not homogeneous, and thus LDA is not exact. For jellium spheres containing only a finite number of electrons, quantum Monte Carlo techniques have been employed [49][50][51] . For infinitely extended jellium on the other hand, going beyond DFT is much more involved 52 and has never been used for the calculation of electronic friction coefficients. Earlier work at the generalized gradient approximation (GGA) to DFT has been done in the context of the effective medium theory (EMT) [53][54][55] .
EMT parameters have been obtained from immersion energies calculated with the atom-injellium model [54][55][56] , and both are modified when using GGA instead of LDA 56,57 . However, the effect of employing GGA instead of LDA on the friction coefficients has not been investigated before. Therefore, friction coefficients are calculated at the GGA level in this work and compared with those obtained with LDA.
Furthermore, spin polarization can affect the value of the friction coefficient for atoms in jellium at low jellium densities, though it is still a matter of discussion whether spin polarization effects should be included within the LDFA scheme when it is applied to molecules 10,56-58 . For carbon it was found that spin-polarized calculations could result in a 70% reduction of the friction coefficient at low jellium density. Nevertheless, this did not alter results for the dissociative chemisorption of methane on Ni(111) 10 . Moreover, a large amount of other elements across the periodic table were found to exhibit a spin moment when spin polarization was allowed within the atom-in-jellium model 57,58 . Furthermore, the jellium can also be spin polarized, to reflect the magnetic moment in ferromagnetic metals 59 and spin friction has also been observed in STM experiments 60,61 . A thorough study into the effect of spin polarization on the friction coefficient will be presented in this paper.
Finally, the atom-in-jellium model is not only important for gas-surface reactions, but also for other kinds of experiments, e.g. analysis of the energy loss of swift (heavy) ions in solids and surfaces [62][63][64][65][66][67][68][69][70][71][72][73][74][75][76] . However, the tabulated data of Puska and Nieminen 45 , that is commonly used in this context, is limited to the first three, incomplete, rows of the periodic table -which is insufficient for studies involving energy dissipation of heavier atoms on metal surfaces [77][78][79] . To extend the amount of elements for which the atom-in-jellium model can be applied, we present here the electronic friction coefficients from hydrogen up to uranium for a variety of jellium densities. Although it is well known that relativistic effects can influence the electronic structure of heavier free atoms 80,81 , to the best of our knowledge friction coefficients have not been obtained whilst employing relativistic LDA. Therefore, we will also investigate the role of relativistic effects for friction coefficients.
The organization of the present paper is as follows: In Sec. II, first the theory behind the atom-in-jellium model is summarized (Sec. II A 1) before relativistic extensions (Sec. II A 2) and computational details (Sec. II B) specific to this paper are described. In Sec. III A, a comparison between the results for electronic friction coefficients obtained with LDA and GGA is made. Section III B concerns spin polarization. Relativistic effects are discussed in Sec. III C. Finally, in Sec. IV we summarize the main conclusions of this paper.

Non-relativistic atom in jellium
The homogeneous electron gas (jellium) is a model for simple metals that consists of a constant positive background and negative electron charge density resulting in an overall neutral system. Both densities are characterized by the density parameter n 0 ≥ 0 a −3 0 and commonly quantified by the Wigner-Seitz radius r −3 s = 4 3 πn 0 , which is the sphere radius of the mean volume of an electron.
Using spherical coordinates, the radial parts of the corresponding continuum of states are given by spherical Bessel functions j l (kr). The (integer) quantum number l ≥ 0 characterizes the angular momentum, whereas the continuous quantum number k ∈ [0; k F ] describes the momentum of the state. The highest occupied state is given by the Fermi energy E F and the concomitant Fermi momentum k F : Summing over momenta and (an infinite amount of) angular momenta yields the electron probability density of jellium which is constant due to l (2l + 1)j 2 l (kr) = 1. Spin-polarized jellium is a simple model for ferromagnetic metals 82 , which introduces homogeneous electron probability densities n σ J (r), σ ∈ {↑, ↓}, in the case of collinear spin considered here, such that The spin-dependent Fermi momenta are given by The strength of the magnetism is characterized by a homogeneous spin polarization ζ = n ↑ J −n ↓ J n 0 , where ζ = 0 corresponds to the original, non-spin-polarized jellium (n ↑ J = n ↓ J = n 0 2 ) and ζ = 1 to the ferromagnetic case (n ↑ J = n 0 , n ↓ J = 0). Throughout the rest of this paper, the spin up channel represents the majority spin channel, i.e., ζ ≥ 0, n ↓ J (r) ≤ n 0 ≤ n ↑ J (r) and k ↓ F ≤ 3 √ 3π 2 n 0 ≤ k ↑ F . In the atom-in-jellium model, homogeniety is destroyed by immersing an atom in a jellium background with density r s . This model can be solved approximately using DFT. Assuming spherical symmetry, the following one-electron Kohn-Sham equations for the radial part of the atom, which is centered at the origin, need to be solved numerically where ψ σ (r) and σ are the radial part and the corresponding eigenenergy for the Kohn-Sham orbitals. Due to the spherical symmetry these orbitals are (2l + 1) degenerate in the (omitted) magnetic quantum number m. The spectrum consists of (localized) bound states ψ b,σ n,l (r), which are characterized by the main (n) and angular (l) quantum numbers, and (delocalized) scattering states ψ sc,σ l (r; k), which yield the total electron probability density analogously to Eq. (3), where n ↑ AIJ (r) = n ↓ AIJ (r) in the non-spin-polarized case. The potential V σ (r) in Eq. (5) is given by where Z is the nuclear charge of the immersed atomic impurity and V σ xc is the exchangecorrelation potential. Choosing V σ xc of the jellium background as the zero reference of the potential [as done in Eq. (7)] yields energy eigenvalues b,σ n,l < 0 (0 < sc,σ l (k) < (k σ F ) 2 ) for the bound (scattered) states. Since V σ depends on the electron distribution, Eq. (5) needs to be solved self-consistently.
The scattering states are normalized by matching them at the cutoff radius R to their asymptotic limit 83 , where j l and n l are the spherical Bessel and Neumann functions, respectively. The phase shift δ σ l (k) is given by where (ln ψ sc,σ The electronic friction coefficient η can be calculated from the difference between the phase shifts δ l (k F ) of the scattering states at the Fermi energy 44,84 : If the jellium background is not spin polarized and the atomic impurity does not induce spin polarization, the summation over the σ in Eq. (11) simply yields a factor two, since the phase shifts for spin up and spin down are identical. Due to the complete screening of the nuclear charge Z by the jellium background, the phase shifts obey the Friedel sum rule 85 , with Z b being the amount of bound electrons. The atom-induced density of states per unit momentum is given by 2. Full and scalar relativistic extension a. RLDA. We have extended the atom-in-jellium model to account for relativistic effects. In the fully-relativistic case the following Kohn-Sham-Dirac radial equations need to be solved 86 , where The zero of the energy is chosen such that = 0 describes electrons with zero kinetic energy in the jellium background (i.e., the rest mass of the electron, c 2 in present units, has been taken out). g(r) and f (r) are the radial parts of the large and small components of the two-component Pauli spinors that describe the Kohn-Sham states, respectively. They are characterized by the relativistic quantum number κ, that is related to the total angular momentum quantum number j = l ± 1 2 according to The potential V R in Eqs. (14) and (15) has the same form as in Eq. (7). In the relativistic local-density approximation (RLDA) used in this paper, a relativistic correction to the (nonrelativistic LDA) is included in the exchange-correlation potential 87 . Again, a self-consistent solution is required because V R depends on the total electron probability density, which is obtained like in the non-relativistic case [see Eq. (6)] as a sum over bound and scattering states resulting from Eqs. (14). For the latter, the boundary conditions of the radial parts of the large and small components are 88,89 respectively, wherel = l − sgn(κ) and The phase shift is then 86,88,89 where This yields electronic friction coefficient η RLDA according to Eq. (11) by summing over κ instead of l and σ.
b. ScRLDA. In addition to the fully relativistic treatment, we have also implemented a scalar-relativistic description according to the approximation proposed by Koelling and Harmon 90 : Eliminating the small component and averaging over the spin-orbit components in Eqs. (14) leads to M ScR is defined analogously to Eq. (15) using the potential V ScR . In the scalar relativistic local-density approximation (ScRLDA), V ScR corresponds to V R , but is based on the electron distribution that is obtained self-consistently with Eq. (21). The total electron probability density is calculated as before [see Eq. (6)] as a sum over bound and scattering states, which are characterized by the same quantum numbers as in the non-relativistic case. For c → ∞ (and thus M σ ScR → 1), Eq. (21) reduces to the non-relativistic case given by Eq. (5). After substituting the corresponding non-relativistic quantum numbers into Eq. (17a), the boundary conditions for the scattering states are identical to the non-relativistic case given by Eq. (8). Consequently, the phase shift is obtained in the same way as in Eq. (9), where the logarithmic derivative (lng sc,σ l ) (R; k) is defined analogously to Eq. (10). The corresponding electronic friction coefficients η ScRLDA can then be calculated according to

III. RESULTS
A. Generalized gradient approximation Figure 1 compares the friction coefficients for Z = 1 − 92 at r s = 2 using LDA and GGA. We reproduce the results presented by Puska and Nieminen 45 for Z = 1 − 18 at various densities and for Z = 1 − 40 at r s = 2 using LDA. The differences between friction coefficients obtained with LDA and GGA are negligible. This is also observed at other jellium densities. The lack of difference between LDA and GGA is also found for the induced density of states. Since the difference between the induced density of states obtained with LDA and GGA is negligible, it is not surprising that the friction coefficients remain unchanged.
This is at odds with the fact that previously it has been reported that including the gradient has an influence on the EMT parameters 56 , specifically the neutral sphere radius and cohesive function, within the same atom-in-jellium model. Puska  On closer inspection, the main correction of GGA over LDA comes from spatial regions where the reduced density gradient ( ∇n(r) n(r) 4/3 ) is large. This correction is particularly relevant wherever the total electron probability density n is low and its gradient is large -as is the case in the exponential tail of the free atom electron density at large distances. Consequently, the exchange-correlation energy and thus the total energy of the free atom is significantly different. Since the latter enters the expression of the immersion energy [Eq. (A1)], GGA yields significantly different values for this EMT parameter.
Friction coefficients on the other hand are entirely defined by the potential that enters the Kohn-Sham equations for the atom in jellium [Eq. (7)].   The differences here are small, ranging from a 15% reduction to 15% increase of the friction coefficients. However, when the background density is lower, spin polarization becomes increasingly more important. Not only are more elements affected by spin polarization, but the differences are relatively larger at lower densities, ranging from a 90% reduction to a 30% increase of the friction coefficients at r s = 3.5. Free atoms with a half-filled d or f orbital are the most affected by spin polarization. At even lower densities (r s > 5) this effect is also observed for half-filled p orbitals. In general, spin-polarized friction coefficients tend to be lower than non-spin-polarized ones. However, a higher friction coefficient is also possible, seen most prominently for free atoms with an almost empty or completely filled orbital.
To understand what is causing the difference between the friction coefficients, we first compare the trends in total spin and the difference between the friction coefficients due to Results obtained with LDA and LSDA are indicated by the red circles and blue crosses, respectively.
Lines are merely to guide the eye. The numerical data is tabulated in the Supplemental Material 95 .
spin polarization across the periodic table in Fig. 4 at r s = 3.5. The appearance of a (nonzero) total spin coincides with the change in the friction coefficient and is only observed at this density for free atoms with partially filled d and f orbitals. The total spin is caused by a difference in the amount of scattering spin-up and -down electrons. The maximum total spin found for atoms with a partially filled f orbital is 3.5. Moreover, for atoms with a partially filled d orbital the maximum total spin is 2.5. At lower density (r s > 5), a total spin for atoms with a p orbital is also observed, with 1.5 being its maximum value. The maximum total spin that is observed in the scattering states corresponds to half-filled f , d, and p orbitals, respectively.
The appearance of a total spin and its effect on the friction coefficient can be understood by looking at the induced density of states [see Eqs. (11) and (13)   coefficient. Once more, this can be understood from the induced density of states. For example, the induced density of states resonance peaks of Cobalt are at a lower energy compared to vanadium. When the resonance peak of Cobalt is split due to spin polarization, the spin down resonance peak is still within the band since the non-spin-polarized resonance peak for Cobalt is at a significantly lower energy than, e.g., for vanadium. The spin down resonance peak, being closer to the Fermi energy than the non-spin-polarized resonance peak, causes a higher induced density of states at the Fermi energy (increase of 60%) and concomitant larger friction coefficient (increase of 30%). The split in the resonance peak near the Fermi energy is also observed for other elements for which spin polarization yields a total spin. Which specific scattering states contribute significantly to the induced density of states close to the Fermi energy, varies with elements and densities. Furthermore, using GGA instead of LDA does not produce different results.
The differences in the friction coefficients using a spin-polarized jellium compared to a non-spin-polarized jellium (ζ = 0) are presented in Fig. 6. As the spin polarization becomes larger, the differences increase as well. Whether the friction coefficient increases or reduces is dependent on the element and density, and as such no clear trend is observed. Again, these differences in the friction coefficients are not caused by the bound states, but by the scattering states. This can also be seen in Fig. 7 where  respect to the non-relativistic friction coefficients of 20% at r s = 1.5. At lower density, these effects are relatively larger, especially for atoms with partially filled d and f orbitals, ranging from a 40% reduction to 180% increase of the friction coefficient at r s = 5. The common trend is that relativistic effects lower the friction coefficient for atoms with partially filled s and p orbitals and increase the friction coefficient for partially filled d and f orbitals.
Another look at the induced density of states is required in order to explain the differences caused by relativistic effects. Figure 10 shows the induced density of states for tungsten (Z = 74), for which relativistic effects increase the friction coefficient and radon (Z = 86), which is affected in the opposite way. In general, the induced density of states at low energies is higher due to relativistic effects. Furthermore, the resonance peak near the Fermi energy is lower and is shifted to a higher energy compared to LDA. How this affects the friction coefficient depends on the induced density of states at the Fermi energy. Typically, the induced density of states will be lower within the ScRLDA if the peak is relatively close to the Fermi energy due to the smaller resonance peak, resulting in a reduced friction coefficient. Otherwise, when the resonance peak is at a comparatively lower energy, the shift of the resonance peak increases the induced density of states at the Fermi energy and the friction coefficient.  is the 5d elements, for which the total spin is partially due to the presence of more bound spin up than spin down electrons originating from the 6s orbital, but this results only in a slight increase of the friction coefficients (< 10%). This effect was not observed with the LDA.
Additionally, in Table I friction coefficients are given for a few heavy elements obtained with ScRLDA and RLDA at r s = 1.5 and 5. Fully relativistic calculations did not alter results significantly compared to ScRLDA. At high density (r s = 1.5), the differences were smaller than 5%. Spin-orbit coupling has a slightly bigger effect (< 10%) at low densities (r s = 5), but the absolute differences at low densities are small, especially compared to the differences between LDA and ScRLDA.  The immersion energy, which describes the energy cost or gain of placing an atom in jellium, is obtained from the atom-in-jellium model by taking the energy difference between the atom in jellium, and the pure jellium and free atom 53,83 , where the energy difference between the atom in jellium and pure jellium can be obtained from a single calculation of the atom in jellium: The difference in kinetic energy is The difference in Coulomb energy is given by and the exchange-correlation energy difference is to calculate immersion energy curves for various first and second row atoms without spin polarization.
Important parameters for the EMT can be obtained from the atom in jellium model 55 .
The so-called cohesive function, Minimizing E c with respect to n 0 yields the cohesive energy E coh = |E c (n coh 0 )| and the concomitant density parameter n coh 0 , which are two very important EMT parameters. We have implemented the calculation of E c and s into LDFAtom, but have not made use of it in the scope of this paper. In our implementation LDFAtom, the radial Kohn-Sham equations [Eqs. (5), (14) and (21)] are solved by rewriting them for the non-relativistic (Schrödinger), RLDA and ScRLDA in the form of two coupled first-order differential equations that are completely equivalent to the respective formulation in Sec. II.

Grids
Using the fourth-order Adams-Bashforth integration method 97 already implemented in dftatom 91 , the equations presented in the preceding Appendix B 1 are solved on a real-space grid, with i ∈ {0, 1, 2, . . . , N } and where is based on the logistic function: This grid enables adequate sampling near the atomic impurity at the origin because the grid points being logarithmically distributed for r 0 ≤ r i < r is . For r is < r i ≤ r N , grid points become more and more equidistant, which adequately samples the long-range part at large distances from the impurity where perturbation of the jellium has (almost) decayed. We have found empirically by extensive convergence tests that α = 36a −1 0 , i s = 2 5 N , and N = 6000 provide a very accurate solution of all calculated properties. After introducing analytic continuations of the spherical Bessel functions j l (kr) for small arguments (kr ≤ 10 −7 ), we have set r 0 = 10 −7 a 0 . r N = R has been varied individually for each atom in a range from 18a 0 to 28a 0 until the Friedel sum rule Eq. (12) is numerically fulfilled within 10 −4 in each case.
A sufficient number of angular momenta (l max ) needs to be included in the calculation of the scattering states, which is ensured by mandating n J (r 0 ) n J (r N ) − 1 < 10 −6 in a separate calculation for the unperturbed jellium background [see Eq. (2)]. Integrations over k [like e.g. in Eq. (6)] are performed with an equidistant grid of 250 points.

Self-consistent solution
For the initial guess of the atom in jellium, the self-consistent density of the free atom is added to the background density of the jellium. The mixing between self-consistent field (SCF) cycles is performed with a limited memory version of Broyden's second method 98-100 .
The self-consistency is evaluated by checking the convergence of the Kohn-Sham effective potential the concomitant eigenenergies. For the former, the Euclidian norm of each spin component of the potential [see Eq. (7)] is calculated at the beginning of each SCF cycle. Likewise, after the potential has been updated to V end (r) at the end of each SCF cycle, the Euclidean norm of the difference with respect to V begin (r) is calculated. If V σ end (r)−V σ begin (r) 2 V σ begin (r) 2 < 10 −6 for both spin channels, then the potential is considered to be (sufficiently) self-consistent. For the Kohn-Sham eigenenergies, only the largest difference between the current and previous SCF cycle is considered and only when the potential already fulfills the aforementioned self-consistency criterion. When this difference is smaller than 5 * 10 −6 a 0 ·Ha, the eigenenergies are considered to be self-consistent as well and convergence is achieved, i.e., the ground-state solution is obtained.
Weakly bound states can cause calculations not to reach self-consistency. This is caused by the appearance and subsequent disappearance of bound states into the continuum between SCF cycles due to the close proximity of these states to the bottom of the continuum at 0 Ha (for the energy zero chosen in LDFAtom, see Sec. II A 1). To stabilize the SCF convergence -i.e., for purely numerical convenience and without any physical meaninga broadening scheme is introduced for the occupation of such weakly bound states using a Fermi-Dirac distribution: f FD ( nl ) = 2l + 1 exp( nl / B ) + 1 . (B8) Here f FD ( nl ) is the occupation number of the bound Kohn-Sham state with energy σ n,l and B is the broadening parameter. The bound state search is stopped when σ n,l > 5 B . We have used 10 −3 Ha < B < 10 −2 Ha and confirmed that this does not affect the friction coefficients significantly. However, even with this approach, SCF convergence could not be achieved in some cases, mainly d and f elements at low jellium densities.