Multiple spin-phonon relaxation pathways in a Kramer single-ion magnet.

We present a first-principles investigation of spin-phonon relaxation in a molecular crystal of Co2+ single-ion magnets. Our study combines electronic structure calculations with machine-learning force fields and unravels the nature of both the Orbach and the Raman relaxation channels in terms of atomistic processes. We find that although both mechanisms are mediated by the excited spin states, the low temperature spin dynamics is dominated by phonons in the THz energy range, which partially suppress the benefit of having a large magnetic anisotropy. This study also determines the importance of intra-molecular motions for both the relaxation mechanisms and paves the way to the rational design of a new generation of single-ion magnets with tailored spin-phonon coupling.


INTRODUCTION
Storing information at the molecular level is a longstanding dream of chemists and material scientists [1]. In magnetic information storage, single-ion magnets (SIMs) based on Kramer ions [2] have received extensive attention due to the possibility to synthesize coordination compounds with large easy-axis magnetic anisotropy and almost-perfect axial symmetry [3,4]. Under these conditions the molecular magnetic moment is stabilized against temperature, since the spin relaxation must occur through the interaction with the lattice vibrations and requires an energy exchange proportional to the size of the molecular magnetic anisotropy, U eff . This relaxation regime is named after Orbach. The relaxation time generally follows an Arrhenius-like law with temperature, τ = τ 0 exp(U eff /k B T ), where k B is the Boltzman constant, τ 0 is the attempt time and T the temperature.
The design of the crystal field of SIMs has been, for decades, the driving strategy to increase U eff [5][6][7][8][9]. Such strategy has recently reached an important milestone with the synthesis of magnetic molecules having a blocking temperature above the nitrogen's boiling point [10]. Despite these accomplishments, the improvement of the spin lifetime of these complexes has not completely matched the expectations raised by designing magnetic anisotropy barriers above 1000 K. This is mostly due to the fact that i) the attempt time τ 0 takes very small values [11,12] and ii) an additional relaxation mechanism, often interpreted as Raman relaxation, appears at low temperature [9]. The two effects pose strong limitations to the progress in this field and, despite their crucial importance, they are still not completely understood at the atomistic level.
Here we provide a microscopic explanation of these phenomena by delivering a comprehensive first-principles description of the spin relaxation of a prototypical SIM with large axial magnetic anisotropy, namely [Co(L) 2 ] 2− , where H 2 L=1,2-bis(methanesulfonamido)benzene [13]. * lunghia@tcd.ie Our method is based on state-of-the-art electronic structure theory [14], first-principles spin relaxation [15] and machine-learning tools [16] and it makes it possible to deliver the a parameter-free description of both oneand two-phonon processes leading to spin relaxation in SIMs, with the inclusion of all the phonons in the Brillouin zone.

RESULTS AND DISCUSSION
The spin properties of transition-metal-based SIMs, such as our Co 2+ complex, are described by a spin Hamiltonian including a single-ion anisotropy term where S describes the |S| = 3/2 ground state of a highspin Co 2+ . The calculated anisotropy tensor, D, correctly captures the easy-axis nature of the complex and predicts an excited Kramer doublet at ∼ 200 cm −1 . Fig. 1A offers a graphical representation of the spin states. In comparison, experiments place the Kramer excited state at 230 cm −1 above the ground state [13]. The spin Hamiltonian of Eq. (1) is the starting point for the treatment of the spin relaxation. The spin-phonon coupling interaction is extracted from the dependence of the anisotropy tensor D on the atomic positions [17]. At the first order in the coupling strength this interaction is described by the Hamiltonian whereQ αq represents the α-phonon with reciprocal vector q and vibrational frequency ω αq . Within this weak-coupling regime, here we consider time-dependent processes at both the first and second order of timedependent perturbation theory.
At the first-order the transition rate between two spin states |a and |b is provided by the expression where G 1−ph = δ(ω − ω αq )n αq + δ(ω + ω αq )(n αq + 1) andn αq = [exp( ω α /k B T ) − 1] −1 is the Bose-Einstein thermal population. The first (second) term of G 1−ph accounts for a spin transition due to the absorption (emission) of a single phonon. For perfectly axial Kramer systems in the absence of external interactions able to break time-reversal symmetry, the direct transition S z = 3/2 → S z = −3/2 is forbidden and relaxation based on Eq. (3) must proceed through the excited state by absorption of a phonon. Fig. 1B schematically reports such a process. At the second order two phonons can simultaneously mediate the spin relaxation, if an intermediate spin state, |c , allows the transition. In this case the rate of population transfer between two spin states |a and |b is given by where G 2−ph ± , similarly to G 1−ph for single phonon processes, accounts for the energy conservation and the phonon thermal population of two-phonon interactions. For instance, a simultaneous absorption and emission of two phonons with virtually the same energy, schematically described in Fig. 1C, would contribute to the T dependence of Here we implement Eqs. (3) and (5) in a first-principles computational scheme able to provide an estimate of the spin relaxation time. The numerical calculation of the spin relaxation time, either at the first or second order of perturbation theory, requires the calculation of the single-ion tensor derivatives, (∂D/∂Q αq ), of the phonons's frequencies, ω αq , and of the normal modes' compositions,Q αq , across the entire Brillouin zone. The numerical calculation of the first-order derivatives of the molecular anisotropy can be performed relatively easily within a first-principles scheme, as discussed before [17], and it requires ∼ 1000 CASSCF simulations. Phonons calculations represent a bigger challenge. SIMs generally comprise tens of atoms and often crystallize together with counter ions in a quite large unit cell. This is also the case for our [CoL 2 ] 2− , whose unit cell contains 396 atoms.
The ab-initio calculation of the phonons at a generic point in the Brillouin zone requires the use of multiples of the unit cell and, except for simple systems [18], it is generally impractical. Here we solve this challenge by employing a machine-learning force field approach. The model is based on the spectral neighbour analysis potential (SNAP) [16] augmented with van der Walls and electrostatic interactions. The latter are based on a Lennard Jones potential and atomic point charges, respectively. After the successful training of the model on just ∼ 200 DFT reference calculations, the force field is used to optimize a 2 × 2 × 2 super-cell containing 3168 atoms and calculate all the lattice harmonic force constants. The machine-learning force field is able to reproduce atomic forces within 3.0 kcal/mol/Å from the DFT reference value. In comparison the training set contains atomic forces that span a 1200 kcal/mol/Å range. A Γ-point DFT calculation of phonons, requiring 2375 DFT runs, is also performed in order to validate the ML model. An absolute mean error of just ∼ 6 cm −1 is observed for vibrations in the energy window of interest. More details on the model validation are provided in the SI. It should also be noted that the training of the ML force field requires a fraction of the computational cost needed to evaluate the sole Γ-point phonons with DFT and massively undercuts the cost of computing phonon dispersions by DFT.
The lattice force constants and the spin-phonon coupling coefficients, (∂D/∂Q αq ), are used together with Eqs. (3) and (5) to compute the spin-relaxation time as discussed in the SI. Results for both first-and secondorder perturbation theory are reported in the top panel of Fig. 2, together with the best fit to experiments [13]. The first-order theory produces the expected exponential dependence of the spin relaxation with T . In contrast, Raman relaxation shows two different regimes: one pseudo-exponential that dominates τ below 30 K and a second one that almost identically mimics Orbach relaxation at high T .
We now proceed to discuss the potential role of lattice anharmonic interactions. Following the basic concepts of reference [19], a finite phonon lifetime can be accounted for by substituting the Dirac functions appearing in G 1−ph and G 2−ph with a Lorentzian function. The predicted spin lifetimes obtained with a constant linewidth ∆ in the range 0.1-1.0 cm −1 are reported in the bottom panel of Fig. 2. The introduction of a finite phonon lifetime only affects the low-temperature regime of the Orbach relaxation mechanism, where τ assumes a T profile virtually identical to the Raman one, with the only difference of being linearly dependent on the parameter ∆. The results for the low T regime, regardless the mechanism considered, are consistent with the commonly observed exponential behaviour, where a polynomial expression is used to model the relaxation time as function of temperature, i.e. τ ∝ T −n , with n often in the range 2 < n < 4 [7,12]. Although the computa-tional results nicely reproduce the experimental profile, the overall τ is underestimated by ∼ 2 orders of magnitude. This discrepancy is attributed to numerical precision. Note that a careful benchmarking and an improvement of the computational methods for obtainingV αq and the phonon-related quantities represent our future methodological challenge.
Let us now discuss the τ (T ) profile in terms of atomistic processes. The slope of ln(τ ) vs T −1 in the Orbach regime, for T > 30K, provides a U eff ∼ 200 cm −1 . This corresponds to the absorption of phonons with energy in resonance to the excited Kramer spin doublet. For a perfectly harmonic phonon bath, this relation is preserved at any temperature. However, the inclusion of anharmonic interactions makes the absorption of outof-resonance low-energy phonons the favourable singlephonon process at low temperature (T < 30K), where the in-resonance phonons are too scarce [19]. Considering the Raman process, the transfer of population within the ground Kramer doublet state is driven by the simultaneous absorption and emission of two phonons with virtually the same energy, as depicted in Fig. 1C. This interaction is made possible by the presence of spin excited states that mediate the interaction. From the dominator of Eq. (5), it is clear that the phonons do not need to be in resonance with the spin excited states. For T < 30K, phonons from the low energy window ( ω αq = 15 − 50 cm −1 ) provide the main contribution to spin relaxation and determine the τ (T ) profile through n αq (n αq + 1) ∼n αq .
The deviation from a perfect exponential regime and the appearance of the polynomial behaviour is understood by considering that by increasing temperature more phonons become populated and start contributing to the spin relaxation. For T > 30K, phonons with similar energy to the excited spin states becomes sufficiently populated and a two-phonon relaxation through the excited Kramer doublet becomes the dominant pathway. This process is initiated by the absorption of a phonon with energy just above the energy of the excited spin states and the emission of a second one that releases the excess of energy. This interpretation of Raman spin relaxation is also consistent with a recent experimental and theoretical investigation of a Lanthanide SIM [20]. It is interesting to note that Raman relaxation, receiving contributions from phonons across the entire Brillouin zone, is always able to fulfil the in-resonance condition with the spin states. For this reason it is not influenced by the inclusion of anharmonic interactions that introduce out-resonance absorption/emission.
An important result of this study is that low-energy phonons, regardless of whether contributing through the anharmonic Orbach or the Raman mechanism, are always able to induce spin relaxation and limit the spin lifetime even at very low temperatures. In this regime, relaxation time still depends on the energy of the excited states, but through the relation τ ∼ 4|D zz − 1 2 (D xx + D yy )| 2 instead of exponentially, as it happens for the high-temperature regime. As a consequence, a strategy solely based on increasing the single-ion magnetic anisotropy might be insufficient to protct spins from these relaxation processes. Besides the coupling of different magnetic centres [21] the engineering of molecular vibrations represent an interesting possibility [22]. First-principles simulations of spin relaxation offer a new perspective on this challenge as they open up the possibility to directly address the origin of the overall relaxation rate, as represented by the spin phonon coupling coefficients (∂D/∂Q αq ). In order to gain an insight on this aspect of the spin dynamics we must unveil the molecular nature of the vibrations contributing to the spin relaxation.
We have then analysed the effective spin-phonon coupling along the direction Γ−X in the Brillouin zone. This quantity is defined as |V αq | 2 = s,t=x,y,z (∂D st /∂Q αq ) 2 . Panels A and B of Fig. 3 report the result of this analysis in the most relevant energy windows for both Orbach and Raman relaxation. Firstly, we note that |V αq | 2 in the high-energy window is one order of magnitude larger that in the low-energy one. This is in agreement with other computational studies, where high-energy optical phonons have been found more strongly coupled to the spin with respect to the acoustic ones and those belonging to low-energy optical branches [15,17]. Indeed, the largest spin-phonon coupling is observed for a mode occurring at ∼ 160 cm −1 . As depicted in Fig. 3D, the atomic displacements associated with this vibration largely involves the modulation of the Co-N bonddistances and, to a minor degree, of all the degrees of freedom associated to the metal-ligand coordination sphere. The strong spin-phonon coupling of this type of motion has also been reported in both theoretical and experimental recent investigations on Co 2+ SIMs [13,23,24]. The vibrational mode that more closely matches the spin excited state energy is pictured in Fig. 3C and it mostly involves the relative rotation of the two ligands with respect to the axis passing through them and the Cobalt ion. Interestingly the first vibrational state at the Γ point corresponds to a rigid molecular translation and it is therefore inactive (see Fig.3E). In contrast, acoustic phonons, generally consisting of rigid molecular translations, get activated towards the Brillouin zone boundary as they couples to molecular rotations [15,18] (see Fig. 3F).

CONCLUSIONS
In conclusion, we have delivered a comprehensive understanding of the spin-phonon relaxation at the atomistic level in a prototypical Co 2+ complex with large anisotropy barrier. We have then provided an interpretation of both the Orbach and Raman spin-relaxation mechanisms. These results suggest that the inclusion of synthetic guidelines that explicitly address the interaction between spins and lattice vibrations have the potential to further push the limits of the spin lifetime in SIMs. In this regards, first-principles spin dynamics sim-ulations provide a unique opportunity to obtain fundamental insights on the spin lifetime in magnetic molecules and drive its optimization.