Zero-point magnetic exchange interactions

Juba Bouaziz ,1,2,* Julen Ibañez-Azpiroz ,3 Filipe S. M. Guimarães ,1 and Samir Lounis 1,4,† 1Peter Grünberg Institut and Institute for Advanced Simulation, Forschungszentrum Jülich and JARA, Jülich 52425, Germany 2Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom 3Centro de Física de Materiales, Universidad del País Vasco/Euskal Herriko Unibertsitatea, 20018 Donostia, San Sebastián, Spain 4Faculty of Physics, University of Duisburg-Essen, 47053 Duisburg, Germany


I. INTRODUCTION
Matter is constituted by a collection of ions and a surrounding cloud of electrons. These microscopic entities obey quantum mechanical laws that, in addition to thermal fluctuations, involve intrinsic quantum fluctuations, a direct consequence of Heisenberg's uncertainty principle. The presence of fluctuations can alter the collective behavior of particles, modifying the physical properties of matter at the macroscopic level, such as the Curie temperature of magnets [1]. In addition, quantum fluctuations determine the energy of the system at its lowest level, the so-called zero-point energy that provides an extra contribution absent in the classical world. A notorious signature is the long-known Casimir effect [2], in which an attractive force emerges spontaneously between two metallic planes separated by vacuum. However, zero-point effects can emerge in a variety of contexts, including recently found light superconducting compounds like LaH 10 , a quantum crystal stabilized by atomic zero-point fluctuations [3], nuclear spin-lattice relaxation of molecular * j.bouaziz@fz-juelich.de † s.lounis@fz-juelich.de magnets [4], and even the internal degrees of freedom in electrical circuit components [5].
As realized in early works [1], quantum fluctuations play a particularly relevant role concerning a central property of the electron, namely, its spin. Known as zero-point spin fluctuations (ZPSFs), they represent an essential ingredient of itinerant electron magnetism and affect a wide range of phenomena, including high-temperature [6][7][8][9] and iron-based [10,11] superconductivity, quantum phase-transitions at 0 K [12], and potentially even skyrmion lattices [13]. Zero-point spin fluctuations are also predicted to play a notorious role in elemental bulk transition-metal paramagnets and ferromagnets [14,15] by modifying the effective magnitude of the spin moments and possibly inducing spin anharmonic effects [16], which can alter the magnetic stability of the ground state [17].
Zero-point spin fluctuations become increasingly important as the size of the system is decreased down to a handful of atoms, a regime where quantum effects prevail. Due to the tremendous current appeal of such low-dimensional systems in the context of information technology and storage as miniaturized magnetic bits [18][19][20], it becomes crucial to understand and control the effect of quantum fluctuations over their magnetic properties. In the ultimate limit of a single magnetic adatom deposited on a nonmagnetic substrate, different measurements display contrasting magnetic trends depending on the probing protocol [18,19,[21][22][23][24]; while ensemble measurements based on x-ray magnetic circular dichroism report the presence of a huge magnetic anisotropy energy (MAE) that protects the magnetic moment against external Distance dependence of the magnetic exchange interaction. (a) Schematic depiction of the magnetic exchange interaction between two fluctuating magnetic moments as a function of interatomic distance (the repeated yellow atom portrays distance dependence). The characteristic RKKY-like dependence is depicted by the oscillating curves, which also exhibit the strong renormalization from the classical (red solid curve) to the quantum (blue dashed curve) prediction. Also shown is the calculated direct exchange interaction as a function of interatomic distance of dimers (b) Co/Pt111, (c) Fe/Cu111, and (d) Fe/Pt111. The bare and fluctuation-renormalized Heisenberg exchange interactions are denoted by red circles and blue triangles, respectively. In (b) and (c) the respective experimental data (black stars) from Refs. [27,28] are also included with permission from Nature Publishing Group.
interactions, local probing techniques generally find no stable magnetic signal in the very same systems. To resolve this apparent contradiction, first-principles theory has successfully invoked ZPSF as a mechanism that destroys the magnetic bistability by locally reducing the MAE barrier; the larger the fluctuations, the larger the local reduction, clarifying many of the observed trends [25,26].
Nevertheless, the magnetic signal of multiatomic quantum devices [19,[27][28][29][30][31] not only depends on the quantum fluctuations of the single magnetic moments, but also emerges from the way magnetic moments talk to each other via the magnetic exchange interaction (MEI). The impact of local and dynamical correlations on the MEI has been investigated by means of dynamical mean-field theory for elemental bulk and surface materials (see, e.g., Refs. [32,33]) without accounting for spin-orbit coupling. However, the impact of ZPSFs on the different components of the MEI in nanoscale systems formed by a handful of atoms remains to be explored, understood, and eventually quantified theoretically and experimentally. Advances in spin-polarized scanning probing techniques undergone in the past decade allowed pioneering magnetometric measurements of the most fundamental MEI between two magnetic adatoms (i.e., a magnetic dimer) [19,27,28]. Remarkably, first-principles calculations based on the local spin-density approximation (LSDA) predict very precisely the nature of the coupling measured experimentally [i.e., ferromagnetic (FM) or antiferromagnetic (AFM)] as a function of the interatomic distance [27,28] [see Fig. 1(a) for a schematic illustration]. However, in direct contrast to this success, the magnitude of the MEI obtained theoretically systematically overestimates the experimental value, with a relative error of more than 100%, as shown in Figs. 1(b) and 1(c).
In this work we show that quantum fluctuations are behind the gap between standard theory and experiments, providing thereby means to quantify them. In particular, we establish a route based on a first-principles dynamical theory for the realistic evaluation of ZPSFs impacting the MEI. By adapting the coupling constant integral formalism [1,7,34] to the modern framework of time-dependent density functional theory (TDDFT) [35][36][37], we show that the reduction of the MEI magnitude induced by local and nonlocal ZPSFs results in very good agreement with previous experimental results [27,28]. In addition, our analysis reveals that antisymmetric spin interactions of Dzyaloshinskii-Moriya type are particularly robust against ZPSFs, implying that quantum fluctuations favor the emergence of chiral magnetic textures. These findings highlight the importance of quantum effects in the study of nanoscale magnets and their future applications. 043357-2

II. ZERO-POINT SPIN FLUCTUATIONS FROM FIRST PRINCIPLES
As the first step in our theoretical analysis, we obtain an expression for the zero-point energy induced by the spin fluctuations, which can then be used to assess the impact on the MEI. For this purpose, we make use of TDDFT linear response theory within the adiabatic LSDA, which represents one of the most powerful tools for analyzing dynamical properties of spins via ab initio methods [35,38]. As a mean-field theory, however, it does not incorporate near-critical fluctuations [39,40], which can be key in low-dimensional systems [25,26]. A notable way in which ZPSFs can be incorporated into the theoretical framework is by making use of the socalled fluctuation-dissipation theorem, a fundamental relation that gives access to the magnitude of the local fluctuation of the magnetic moment [7,25] with χ +− i (ω) representing spatial average of the interacting transverse susceptibility of a magnetic atom at site i [36,37], evaluated using the Dyson-like equation Above, χ +− 0 (ω) and K represent the Kohn-Sham susceptibility, describing electron-hole excitations and the exchangecorrelation kernel, respectively; the underline is used to denote the tensorial character of objects. Going one step further, we obtain the zero-point energy associated with the ZPSFs by combining the modern TDDFT framework with the coupling constant integral method [7,34,41], giving rise to the first important relation (see Appendix A for the detailed derivation) where the trace runs over the number of magnetic units present in the system. We now introduce the fluctuation-corrected magnetic exchange interactions J SF i j , which are defined assuming an extended Heisenberg Hamiltonian of the form Here E SF represents the fluctuation-corrected band energy of the system and E MAE the local contribution from the on-site MAE. In the first-principles framework, E SF is given by the sum of the LSDA band energy E b and E ± of Eq. (3). Then the magnetic exchange interactions renormalized by quantum fluctuations can be obtained from the curvature of E SF as with the bare exchange interaction parameters that are commonly computed using standard DFT. In Eqs. (5) and (6), e α i denotes the transverse component α of the vector e i defining the orientation of the magnetic moment at site i.
Therefore, the challenge consists in finding an expression for the extra term − ∂ 2 E ± ∂e α i ∂e β j in Eq. (5). Relying on the magnetic force theorem together with the infinitesimal rotational method [42][43][44], the correction of the elements of the MEI tensor is determined as (see Appendix C) As a final step, we derive a simple expression for the renormalized MEI in terms of the fluctuating moments by mapping the ab initio χ +− (ω) to the one obtained from the Landau-Lifshitz-Gilbert (LLG) model. This allows the identification of the analytical dependence between the dynamical susceptibility and the MEI, giving rise to (see Appendix C) where ξ 2 i j,± denotes a nonlocal contribution to the ZPSF, Equation (8) represents a central result of the present work, as it provides a quantitative expression showing how the standard tensor of MEI is renormalized by quantum spin fluctuations. As a general and most important trend, this equation shows that fluctuations systematically reduce the predominant symmetric component of the tensor of MEI, given that the hierarchy |ξ i,± | < |M i | holds in general; we will refer to this important contribution (J αα i j ) as the Heisenberg exchange interaction (HEI) component. On closer inspection, Eq. (8) further reveals that the antisymmetric piece of the tensor (J SF,αβ i j with α = β) giving rise to the Dzyaloshinskii-Moriya interaction (DMI) increases upon the action of ZPSF due to its antisymmetric nature, i.e., the property J αβ i j = −J αβ ji . The combination of these features may have profound implications for the stability, wavelength, and shape of intensively studied chiral spin textures, such as spirals, skyrmions, and antiskyrmions [45], which are strongly dependent on the ratio of DMI and HEI.

III. QUANTUM CORRECTIONS AGAINST EXPERIMENTAL EVIDENCE
With the purpose of connecting the developed firstprinciples framework to the microscopic experimental designs of Refs. [27,28], we investigate magnetic dimers deposited on nonmagnetic metallic substrates as a function of interatomic distance, as schematically depicted in Fig. 1. This analysis exposes the oscillatory behavior of the MEI, bringing into play different exchange mechanisms. In the short-range limit, direct exchange dominates due the strong hybridization among the adatom's d orbitals, while for larger distances the 043357-3 Ruderman-Kittel-Kasuya-Yoshida (RKKY) mechanism prevails, indirectly mediating the interaction via the conduction electrons of the metallic substrate [46][47][48].

A. Computational details
We consider Co and Fe dimers deposited on the fcc stacking sites of the Pt(111) and Cu(111) surfaces, respectively; these systems were assessed experimentally, via scanning tunneling microscope (STM) measurements, and theoretically, via bare LSDA calculations, in Refs. [27,28]. In addition, we further investigate an fcc-stacked Fe dimer on Pt(111) in order to provide theoretical predictions on the fluctuation-renormalized MEI that can be tested in future STM experiments. The dimers are relaxed by 14% and 20% for Cu(111) and Pt(111) surfaces [49,50]. Briefly, we first determine the nature of the coupling (i.e., FM or AFM) using the infinitesimal rotation method, while the easy axis (MAE) is resolved from the static transverse susceptibility [51]. Subsequently, the spin-excitation spectrum is computed from the corresponding collinear ground-state configuration, in the spirit of the magnetic force theorem approach.
The simulations have been carried out using the scalarrelativistic all-electron Korringa-Kohn-Rostoker Green'sfunction method, including the spin-orbit interaction selfconsistently [52,53]. An angular momentum cutoff of l max = 3 and a k mesh of 600 × 600 have been used for the construction of the Green's functions in real space. The magnetic impurities have been embedded into a real-space impurity cluster containing 56 sites and 30 substrate atoms. The magnetic excitations have been assessed in the framework of TDDFT [35,37], using the adiabatic LSDA for the exchangecorrelation kernel. Finally, the ZPSFs have been obtained from the frequency integral of the imaginary part of the magnetic susceptibility, with a cutoff frequency set at 250 meV (given that the spin-excitation energies are around few meV), after which an ω −2 decay is assumed [25,26].

B. Results
We begin by briefly describing the calculated groundstate magnetic properties of the adatoms conforming the dimers, which evolve large magnetic moments of 2.3μ B for Co/Pt(111), 3.3μ B for Fe/Cu(111), and 3.5μ B for Fe/Pt(111). In all the systems, the MAE is of the order of a few meV and favors an out-of-plane orientation. These values are found to be mildly affected by the interatomic distance.
Taking Co/Pt(111) as a case study, next we illustrate the main properties of the spin-excitation spectrum of the dimer as a function of interatomic distance in Fig. 2(a). Our calculations reveal that the spectrum is largely dominated by a strong acoustic mode at ∼4 meV. The optical mode (not shown) is much weaker, lying above 60 meV for the shortest distance and quickly merging with the acoustic mode as the atoms are moved far apart. We note that the position of the acoustic mode correlates with the magnitude of the MAE, while the broadening of the spin excitation is induced by electron-hole excitations [35,37]. Our results evidence two clear regimes as a function of interatomic distance: the dimerlike ( 0.7 nm) and the atomiclike ( 0.7 nm). In the dimerlike regime, the symmetry is lowered [51] and the interaction between the two magnetic adatoms makes the acoustic mode highly distance dependent, converging to the single-adatom limit (atomiclike regime) as the strength of the interaction decays.
Owing to the relation established by the fluctuationdissipation theorem [cf. Eq. (1)], these two markedly different regimes are translated into the evolution of the local and nonlocal ZPSF depicted in Fig. 2(b); while the size of the dominant local spin fluctuations oscillates between ∼1.3μ B and ∼1.5μ B at short distances, it converges to a steady value 043357-4 of ∼1.4μ B in the atomiclike regime. In turn, the nonlocal ZPSFs decay quickly as a function of distance from ∼0.7μ B to virtually zero and are therefore only relevant in the dimerlike regime. Notably, the calculated magnitude of the ZPSF is of the same order as the ground-state magnetic moment itself [also depicted in Fig. 2(b)], a huge relative value for a purely quantum effect.
Having determined the evolution of the spin-excitation spectrum and ZPSFs as a function of interatomic distance (Fig. 2), we now assess the impact of the calculated spin fluctuations on the renormalization of the HEI, namely, the average of the diagonal components of J SF,αβ i j in Eq. (8). Calculated results for Co/Pt(111), Fe/Cu(111), and Fe/Pt(111) are presented in Figs. 1(b), 1(c), and 1(d), respectively; for the two first systems, we additionally include experimental STM data available from Refs. [27,28]. The adatoms couple ferromagnetically (positive HEI) for the nearest-neighbor distance and display a characteristic RKKY decaying oscillatory behavior as the interatomic distance is increased, with different oscillation periods determined by the Fermi surface of the substrate. Notably, the renormalization effects caused by quantum fluctuations fix the disagreement between standard theory and experiments at virtually all interatomic distances, as demonstrated in Figs. 1(b) and 1(c). As revealed by our results, the main role of the fluctuations is to systematically reduce the magnitude of the HEI, correcting for the overestimation arising from the LSDA. We note that in Refs. [27,28], the bare HEI was systematically divided by a factor equal to either 3 or 2 in order to theoretical results closer to experiments. Remarkably, including ZPSFs achieves an excellent comparison without the need of invoking a posteriori parameters, thus proving that quantum fluctuations are a fundamental mechanism at play in these type of low-dimensional magnets.

IV. UNDERSTANDING QUANTUM FLUCTUATIONS VIA THE LLG MODEL
In order to identify the role of the fundamental factors that determine the behavior of ZPSFs, we resort to the widely used LLG model for characterizing the microscopic spin dynamics [54]. Our main goal is to obtain an analytical expression of the TDDFT spin susceptibility entering Eq. (7). Then, by working out the evolution of the magnetic moments forming a dimer, we will synthesize the main properties of the local and nonlocal ZPSFs in terms of the physical parameters entering the LLG model.
The LLG equation describing the damped precessional motion of two magnetic moments M i (i = 1, 2) on top of a substrate is with γ the gyromagnetic ratio and η the Gilbert damping. The first term on the right-hand side represents the torque generated by an effective field B eff The second term on the right-hand side of Eq. (10) models the damping process that drives the magnetization back to equilibrium.
The low symmetry of the problem (C s ) dictates that the quantities involved in the LLG equation for a magnetic dimer have tensorial form (see Appendix B). However, for the sake of clarity, in the following we consider a simpler model with M 1 = M 2 ≡ M and where the ZPSFs are determined as a function of two central quantities. The first one consists of an effective MEI weighed by the magnetic moment, J eff = J i j γ /M. We note that this quantity accounts for the distance dependence of the system, given that a large (small) J i j corresponds to a small (large) interatomic distance. The second ingredient is the Gilbert damping η, a quantity that is closely connected to the width of the spin-excitation peak [55] shown in Fig. 2(a) and is known to be a key player for local ZPSFs, as it quantifies the magnitude of electron-hole Stoner excitations [25]. In the dimers studied in this work, the calculated values of η range between ∼0.1 for Fe/Pt(111) and ∼0.4 for Co/Pt(111).
The key information is encoded into the local and nonlocal LLG spin-flip susceptibilities, which can be synthesized as with K eff = Mγ /M, where M denotes the MAE. The ZPSFs are obtained via the frequency integral of the imaginary part of Eq. (11), which we have numerically computed using the trapezoidal rule in terms of the quantities J eff , K eff , and η. Figure 3 illustrates the LLG solutions for the local and nonlocal contributions to the spin-fluctuation amplitude as a function of J eff and η in the region relevant for experiments (we have set the MAE along the out-of-plane direction, matching the experimental situation). The figure reveals valuable information on the nature of the ZPSFs. First and foremost, quantum fluctuations induce a much larger impact on the AFM regime (J eff < 0) as compared to the FM one (J eff > 0). This result holds for both local and nonlocal contributions and supports the generally assumed notion whereby ferromagnets are more robust against quantum fluctuations than their antiferromagnetic counterparts [51,56]. As a second major message, Fig. 3 shows that the magnitudes of both the effective MEI and Gilbert damping play a different role depending on the nature of the magnetic coupling. In the FM regime, η tends to significantly increase both local and nonlocal quantum fluctuations (in line with a previous study on the single-adatom case [25,26]), whereas the effect of J eff is much milder. In turn, the AFM regime shows a more complex pattern. The effect of |J eff | is significant, inducing large local and nonlocal fluctuations as it increases; as for the Gilbert damping, it induces larger local fluctuations but reduces the nonlocal ones, resulting in a competing mechanism. As an important remark we note that for J eff 0, the nonlocal ZPSFs vanish in both the FM and AFM regions independently of η, 043357-5 a feature that explains the quick fall of ξ 2 i j,± observed in the ab initio results of Fig. 2(b).
In order to extract a final lesson, let us focus on the limit of zero damping, i.e., η → 0. This allows the identification of the intrinsic local ZPSFs, which are predominant at large interatomic distances. In this regime, we have ξ 2 ±,i γ M i /2, which is the lower limit of the ZPSFs [see vertical scales in Figs. 3(a) and 3(b)] and provides a similar but inferior approach for the renormalization of the HEI components from Eq. (8): Notably, this simple expression can be readily incorporated into the conventional DFT framework when computing the MEI with virtually no added cost, given that it only involves ground-state magnetic moments. Note also that Eq. (12) provides a sensible macroscopic limit, whereby the renormalization due to quantum fluctuations tends to vanish for large values of the magnetic moment.

V. CONCLUSION
Besides providing insight into the features of the ab initio calculations, the parameter-space map built in Fig. 3 provides a simple recipe to design robust collinear nanomagnets. As its central prediction, an optimal shield against quantum spin fluctuations can be achieved by a combination of large ferromagnetic HEI coupling, relatively low Gilbert damping (i.e., low electronic hybridization), and a strong out-of-plane MAE. On the other hand, in remarkable contrast to the direct-exchange mechanism, ZPSFs play in favor of the antisymmetric nature of the DMI, as proven by Eq. (8) (see also Appendix C). This finding provides an important pre-diction, namely, that quantum spin fluctuations enhance the noncollinearity and the stability of chiral magnetic textures such as skyrmions, governed by the DMI-to-HEI ratio that increases upon the action of ZPSF. This counterintuitive effect appears to be particularly relevant for the ongoing miniaturization process of skyrmions, whose characteristic length scale is approaching the size of the lattice constant and classical predictions start to break down [57].
In summary, we have analyzed the role of quantum spin fluctuations in the fundamental magnetic coupling mechanism of atoms within a TDDFT framework. We have furthermore provided a simple prescription for quantitatively assessing the impact of quantum spin fluctuations in the renormalization of magnetic exchange interactions [cf. Eq. (12)], an expression that can be readily applied to the widespread DFT framework. By applying the approach to magnetic dimers deposited on metallic substrates, we demonstrated that accounting for zeropoint spin fluctuations results in better agreement with the experimentally available data, fixing a relative error of a factor 2-3. Our work therefore suggests that ZPSFs are an important ingredient for the study of nanomagnets and their applications [29,30].

ACKNOWLEDGMENTS
We are grateful to Jens Wiebe and Roland Wiesendanger for sharing with us the experimental data published in Refs. [27,28].  [58].

APPENDIX A: TRANSVERSE SPIN FLUCTUATIONS AND ZERO-POINT ENERGY
In this Appendix we derive the zero-point energy present in the system due to ZPSFs, which emerge as a consequence of the electron-electron interaction. For clarity, a single-band Hubbard model is employed [59], where the Hamiltonian can be written aŝ where the first termĤ 0 accounts for the electrons hopping t i j from site i to j and the second oneĤ int contains the effective Coulomb electron-electron interaction U i . In addition,n iσ = c † iσ c iσ is the occupation operator, with c † iσ (c iσ ) the creation (annihilation) operator of an electron at site i with spin σ . The sum runs over all the sites of the system. Using the commutation relations between c iσ and c † iσ , the interaction HamiltonianĤ int can be reexpressed in terms of the total occupation numbern i =n i↑ +n i↓ and the net magnetic moment m iz =n i↑ −n i↓ as [60] , Eq. (A2) can be recast in terms of spin operators The resulting interaction energy E int can be obtained using the coupling constant integral method [7,34,41]. Using λ as a constant parameter that scales U i adiabatically, The fluctuation energy E is defined [7] as the difference between the total E int calculated in Eq. (A4) and the Hartree-Fock energy (mean-field energy), denoted here by · · · 0 . Therefore, where E includes longitudinal and transverse spin fluctuations, as well as charge fluctuations. We focus our subsequent analysis on the transverse spin-fluctuation contribution alone for two main reasons: First, it represents the largest contribution by far for low-dimensional magnets [25], and second, it is the relevant piece for the analysis of the magnetic stability since its energy scale coincides with that of the magnetic anisotropy energy (meV). Consequently, we focus on the transverse spin-fluctuation energy denoted by E ± and given by This equation can be rewritten using the spin-flip susceptibilities χ +− i (ω) and χ −+ i (ω) via the fluctuation-dissipation theorem [61] The last line has been obtained using the symmetry relation [χ −+ (ω)] * = χ +− (−ω). The transverse spin-fluctuation contribution to E ± [Eq. (A6)] can be recast into [7] E ± = 1 2π Im Tr Above, the trace runs over the number of atomic sites present in the system. In addition, U is a diagonal matrix containing the elements U i and χ +− 0(λ) (ω) is the Hartree-Fock (λU iinteracting state) susceptibility. In the random-phase approximation we have χ +− λ (ω) = χ +− 0 (ω) + χ +− 0 (ω)λU χ +− λ (ω) and E ± finally reads Notably, this expression can be used in the framework of TDDFT by making the connections [36] χ +− 0,i j (ω) → χ +− KS,i j (ω) = dr dr χ +− KS,i j (r, r , ω), with χ +− KS the Kohn-Sham susceptibility and K i the adiabatic exchange-correlation kernel, which is local in space. Note that the radial integrals are performed around each atomic site i or j.

APPENDIX B: MAPPING TDDFT TO THE LANDAU-LIFSHITZ-GILBERT EQUATION
In order to interpret results obtained in the framework of TDDFT calculations, we consider a Heisenberg model for localized magnetic moments, whose dynamics are dictated by the LLG equation [54] Above, the first term accounts for the precession of the magnetic moment while the second one takes into account relaxation effects, with η defining the Gilbert damping parameter. The effective field B i,eff (t ) = B i,ext (t ) + B i (t ) includes the transverse time-dependent external magnetic field driving the magnetization out of its equilibrium orientation B i,ext (t ) in the Kohn-Sham susceptibility χ +− KS (ω) with respect to a modification of the α component of the vector where χ +− (ω) is the spatial average of the spin-flip susceptibility computed within TDDFT and K is assumed to be independent of e α i (i.e., isotropic). Note that the expression (C1) was obtained in the limit of magnetic anisotropy energy K i and magnetic exchange interaction J i j K i , which permits one to disentangle the collective spin excitations and the Stoner contributions to the ZPSF energy. Indeed, K i is on the order of eV, while the energy scale of the spin excitations is the range of a few meV, as studied in the present work. In addition, as we checked analytically and numerically, the frequency decay of the spin susceptibility [61] leads to vanishing ZPSF energy contributions at high frequencies. Thus, the variation of the integrand defined in Eq. (A9) due an infinitesimal rotation of e i is written as To access the magnetic interactions (i.e., curvature of the energy [42]) a second differentiation with respect to another independent variable e β j is performed: The frequency integral of this equation leads to the ZPSF correction of the exchange interactions: Assuming that the electron-hole excitations are weakly affected by the infinitesimal rotation of the magnetic moments, the first term on the right-hand side of the preceding equation is found proportional to Having obtained the general form, we next focus on the particular system analyzed in the main text constituted of two magnetic atoms sitting on a nonmagnetic surface. In this case, the trace in Eq. (C5) becomes where we have introduced χ +− αβ,i j = where ξ 2 i,± designates the aforementioned square of the local spin fluctuation amplitude of the moment M i [25], and the nonlocal contributions ξ 2 i j,± are defined as Note also that the nonlocal spin fluctuation contribution that corrects J αα 12 (for α = x, y) in Eq. (C7) is proportional to J zz 12 . In general, these components may differ due to the anisotropy of the system. However, these differences are relevant in the investigated cases when the atoms are far apart [62], where ξ 2 12,± is vanishingly small (as shown in the main text). In this case, Eq. (C7) simplifies to