The role of anharmonic phonons in under-barrier spin relaxation of single molecule magnets

The use of single molecule magnets in mainstream electronics requires their magnetic moment to be stable over long times. One can achieve such a goal by designing compounds with spin-reversal barriers exceeding room temperature, namely with large uniaxial anisotropies. Such strategy, however, has been defeated by several recent experiments demonstrating under-barrier relaxation at high temperature, a behaviour today unexplained. Here we propose spin–phonon coupling to be responsible for such anomaly. With a combination of electronic structure theory and master equations we show that, in the presence of phonon dissipation, the relevant energy scale for the spin relaxation is given by the lower-lying phonon modes interacting with the local spins. These open a channel for spin reversal at energies lower than that set by the magnetic anisotropy, producing fast under-barrier spin relaxation. Our findings rationalize a significant body of experimental work and suggest a possible strategy for engineering room temperature single molecule magnets.

SUPPLEMENTARY NOTE 1. STOCHASTIC MODELING OF PHONONS LINE-WIDTH.
The diagonal elements of the spin reduced density matrix can be expressed through the second order time convolution less master equation provided the validity of the weak coupling regime for the spin and the phonons degrees of freedom (Born approximation). Under this approximation, Supplementary Eq. 1 can be used to describe the non-markovian dynamics of the spin open-system. The quantity G(t, ω ij , ω α ) is the α-phonon green function calculated at the time t and it is defined as where the subscript 'B' stands for all the degrees of freedom of the phonon bath and the operator T r B ...ρ eq B represents the ensemble equilibrium average over the bath. In Supplementary Eq. (2)q α = 1 √ 2 (â † α +â α ) is the displacement operator for the α-th mode withâ † α (â α ) being the creation (annihilation) operator for that mode of frequency ω α , andρ eq B is the equilibrium density matrix for the bath. G(t, ω, ω α ) contains all the information on the bath dynamics, time and temperature dependence included. This mathematical object is usually studied perturbatively, where the harmonic approximation is used as the ground for the Dyson's equation and the anharmonic interactions between phonons represent a perturbation acting on the normal modes. The presence of anharmonic interactions makes the life-time of phonons finite and consequently the bath spectral shape broadened. In order to access this information about the bath, one requires calculating the anharmonic frequencies and the phonon self-energy [1]. In this work we have tackled the problem through a stochastic approach, which has the advantage of not requiring information about the single-phonon phonon-phonon interactions. Within this framework, the bath Hamiltonian is assumed to have an harmonic form,Ĥ ph = αĤ ph,α = α ω α (n α + 1 2 ), and all the anharmonic features are introduced in a statistical way, according to the fact that the interaction of a mode with the bath (made by all the other modes) is the same interaction that drives it to the thermal equilibrium described by the canonical distribution. If a phonon is at thermal equilibrium, its energy will experience fluctuations and consequently, any eigenstate ofĤ ph will have a finite life-time. Clearly there is an equivalence between the statistical picture and the mechanistic microscopic description of the phonon-bath system [2]. In order to bridge the statistical description with the microscopic one, the phonon Hamiltonian is considered to be time dependent with a time propagator e −i t 0 dτĤ ph (τ ) . The expression for the time-dependent bath normal mode displacement operators iŝ q(t) = e i t 0 dτĤ ph (τ )q e −i t 0 dτĤ ph (τ ) .
Therefore we have The time dependence ofĤ ph enters into ω α that can be expressed as ω α (t) = ω α + δω α (t). According to the Born approximation (ρ B (t) ∼ρ eq B ), we assume that the bath relaxes to equilibrium as its equilibrium fluctuations do, in a fluctuation-dissipation theorem fashion. Therefore, the phonons' energies will be considered randomly fluctuating around the microcanoical (NVE) ω α frequencies with a gaussian distribution of amplitude ∆ = δE 2 = E 2 − E 2 , where all the mean values are taken at the canonical ensamble (NVT) equilibrium. As the time dependent part of the bath Hamiltonian is not due to its operator part we obtain [Ĥ ph (t),Ĥ ph (t )] = 0 and thus e i t 0 dτĤ ph (τ ) = t n=0 e i∆τĤ ph (τn) , where all the exponents commute. According to this, it is then possible to demonstrate that so that Supplementary Eq. (6) becomes Before proceeding further we should consider the relation Finally, according to Kubo's model [3,4], we can write We can then expand F (t ) as a cumulant expansion of averages and take advantage of the Gaussian statistic of fluctuations which, by making the substitution τ = τ − τ , becomes Finally, by assuming that the fluctuation correlation function decays exponentially with a characteristic time τ α i.e. δω(τ )δω(0) = ∆ 2 α exp(−τ /τ α ), we obtain The Heisenberg principle offers a way to estimate the maximum limit of the τ α according to a given energy fluctuation magnitude ∆ α , i.e. τ α = ∆ −1 α . According to its statistical mechanics definition, ∆ α can be evaluated as Supplementary Eq. (15) describes the temperature dependence of the phonon spectra line-width. At low temperature ∆ converges to zero asymptotically and a 1/β linear divergence is observed for β → 0. A similar temperature behaviour is experimentally observed through the determinations of the homogeneous contribution to the phonon line width [5][6][7]. Although Supplementary Eq. (15) qualitatively describes the correct observed behavior for ∆, a few discrepancies between this model and experiments exist. First of all, the experimental line-width approaches a finite value different from 0 for T → 0. Moreover, this model treats on the same footing every modes while experimentally the magnitude of the ∆ divergence depends on the anharmonicity of the mode considered. However, given its extreme simplicity, the proposed model represents an appealing way to introduce anharmonic contributions in the bath Green's function without any further requirements apart from harmonic modes calculation. Finally, it is possible to define the phonon spectral shape according to the phonon Green's function wheren α = 1 e β ωα −1 is the phonon occupation number according to the Bose-Einstein statistics. Supplementary Eq. 17 can be further simplified under the condition that the integrand function would decay on a faster time-scale than the one under which the spin life-time evolves (Markov approximation). Under this assumption, the decay of F(t) becomes purely exponential and the integration limits of Supplementary Eq. 17 can be replaced with ±∞ to obtain a time independent bath Green function G(ω, ω α ): By using ∆ α as defined in Supplementary Eq. (15), the life-time of the normal mode ω=36 cm −1 is in the ns -ps range for temperatures in the 2 − 10 K interval. This is orders of magnitude shorter than the spin relaxation time observed for [(tpa P h )Fe] − (ms -µs), demonstrating that Markov's approximation is valid for the dynamics discussed in the main text of the paper. When the Markov's condition is not fulfilled, one has then to replace Supplementary Eq. (21) with the more general Supplementary Eq. (17).
The three-level S = 1 model has been discussed in the main text and it is defined by the following energy ladder: As discussed in the main text, we assume that there is only one phonon able to directly interact with the spin system. The eigenstates populations p n (n = 0, 1, 2) evolve in time according to the kinetic equations If one assumes an identical spin-phonon coupling strength V , for all the transitions, the various kinetic constants are defined as The kinetic constants k 3 , associated to the transitions between the quasi-degenerate lower-lying states, |0 → |1 and |1 → |0 , have been set to V [2n( ω) + 1]G(δ, ω). This is because of the assumption of E 1 − E 0 = δ ∼ 0 and the fact that the phonon has a frequency ω > 0 and a finite line-width. In this scenario the transition between the two almost degenerate states is induced by the tail of the phonon distribution at an energy ∼ 0.
If we now take the quantisation axis along the direction of an external magnetic field, or equivalently with the z-direction oriented in order to reduce the rhombic terms in the spin Hamiltonian, the magnetisation will be in good approximation M z (t) = p 0 − p 1 . Thus the dynamics equation for M z can be written by subtracting the first and the second equation in (22), namely From Supplementary Eq. (26) it appears clear that M z (t) decays exponentially in time with the decay coefficient determined by the two kinetic constants k 3 and k 1 . In order to understand the nature of the relaxation processes related to these two constants it is necessary to obtain eigenvectors and eigenvalues of the entire system of equations (22). However, instead of solving the Eqs. (22), it is more instructive to separate the problem into the direct intra ground-state doublet population transfer and the excited state mechanism (Orbach's mechanism). This is done, respectively, by forcing k 1 = k 2 = 0 and k 3 = 0. Case 1. The solutions of the system (22) for the k 1 = k 2 = 0 case are and thus Under these conditions the magnetisation vector relaxes to its equilibrium value with a relaxation time τ = 1 2k3 . In the anharmonic case, i.e. assuming ∆ as form Supplementary Eq.15, and in the ω k B T δ limit we have and so the relaxation time follows an Arrhenius like-temperature dependence In the k B T ω case we have In the harmonic case G( ω, δ) = δ( ω − δ) and therefore the relaxation would proceed only when a very-low-energy phonon is available to couple with the spin system. In this situation the most interesting case is k B T δ and under these assumptions the population transfer rate reads where the linear dependency of the direct relaxation mechanism with respect to T is recovered, i.e.
Case 2. When k 3 = 0, a situation describing a pure Orbach's relaxation mechanism, we have two eigenvalues beside the null one: k 1 and 2k 2 + k 1 . The eigenvector corresponding to k 1 is the same as for the case k 1 = k 2 = 0 just discussed and it corresponds to a population transfer from |0 to |1 or vice versa: Interestingly, this process, which introduces the k 1 kinetic constant in the magnetisation decay profile, does not transfer population to the excited states, but involve a virtually instantaneous excitation and emission of the system. It must be noted that even though the population of the excited state |2 is not involved in this process, the presence of this state is indirectly accounted for in the definition of the kinetic transfer rate constant k 1 . The effect on the system of the eigenvalue 2k 2 + k 1 is different from the former one and it corresponds to an eigenvector of the form This kind of eigenvector does not produce relaxation but merely adjust the population of the excited state to its equilibrium value and indeed it does not appear in the M z (t) rate equation (26). Focusing on the sole through excited state mechanism, the transition rate is which, in the resonance condition, becomes for all temperatures. However, in the low-T regime, the general expression (36) implies a relaxation time constant It must be noted that when taking the definition of the phonon life-time of Supplementary Eq. (15) we are implicitly assuming that the magnitude of the phonon frequency sets the temperature at which the anharmonic effects become visible. Due to this feature of the model the phonon damping and its frequency are correlated quantities and it is no longer possible to recover the limit of zero damping without changing the harmonic features of the vibrational mode. For this reason the harmonic limit cannot be recovered from Supplementary Eq. (36) but must be taken before substituting Supplementary Eq. (15) in the definition of k 1 . With this precaution, the population transfer rate in the harmonic limit is As for the direct mechanism, the relaxation in the harmonic regime is possible only under the resonant condition ω = U 0 . In the limit ω k b T the population transfer rate is and therefore where the usual identity U eff = U 0 is recovered.

SUPPLEMENTARY NOTE 3. SPIN HAMILTONIAN PARAMETERS CALCULATION.
Let us consider the system to be described by an electronic HamiltonianĤ =Ĥ BO +Ĥ SO , whereĤ SO refers to the spin-orbit component of the Hamiltonian andĤ BO refers to the Born-Oppenheimer one. The mapping between the electronic Hamiltonian,Ĥ, and the spin Hamiltonian,Ĥ S , is made through the relation where |SM S are the spin eigenstate of total spin S and its z-component M S . The (2S+1) lower-lying eigenstates of H could be written as function of their eigenvalues and eigenvectorŝ Therefore, the mapping with the spin Hamiltonian,Ĥ S , assumes the form Once E k and |k are obtained by means of electronic structure calculations, the system of linear equations (44) could be solved by means of a least square fit [8]. Although this relation is not a priori valid and the quality of the fitting should be evaluated each time, its range of applicability is wider with respect to the perturbation treatment of relativistic interactions [9]. For instance, it can also be exploited when higher order spin Hamiltonian terms are needed to describes the system; this feature is particularly appealing as analytical expressions are missing for these terms. At the same time, a clear connection between specific interactions included inĤ and spin Hamiltonian terms is here lost and the correctness ofĤ S choice falls under the quality of the fit itself. Over parametrisation problems could arise, especially when multi spin complexes are considered. This procedure could be also applied to calculate the spin-phonon coupling coefficients in a complete general fashion. Considering an effective spin HamiltonianĤ S = kq B kqÔkq ( S) [10], expressed as an expansion over even-order tesseral spin operators,Ô kq , theĤ 0 matrix element derivatives could be evaluated as follows: In deriving Supplementary Eq. (45) we take advantage of the fact that the spin part of the wave-function does not depend on the bath normal modes, and of the spin hamiltonian definition b|Ĥ 0 |a = b|Ĥ S |a . The expression (45) can be readily applied to any spin Hamiltonian form.

SUPPLEMENTARY NOTE 4. AB-INITIO CALCULATIONS DETAILS.
Spin dynamics calculations have been carried out for the anion [(tpa Ph )Fe] − SMM [11], where tpa Ph = tris(pyrrolylphenyl-methyl)amine is the ligand of the S = 2 high-spin Fe 2+ ion. The trigonal pyramidal coordination geometry of the amine ligands around the Fe 2+ center makes the first excited S = 2 state almost degenerate with the ground state. This situation makes the spin-orbit interaction especially effective to split the degeneracy of S = 2 ground state multiplet, producing a remarkable easy axis anisotropy with a zero-field splitting of D = −27.5 cm −1 and with an almost vanishing rhombicity. [(tpa Ph )Fe] − crystallizes in a P1 space symmetry group with a primitive cell that comprises two SMMs units and two co-precipitate Na ion coordinated by three dimethoxyethane solvent molecules each, for a total of 228 atoms. The [(tpa Ph )Fe] − structure has been optimized both for the isolated molecule and for the periodic unit cell in the gamma-point approximation. From now on these two models will be denoted as Isol Model and Bulk Model, respectively. The structures produced by the two optimisation schemes are very similar, as expected from the high stiffness of the ligands cage. Supplementary Table 1 reports the first coordination shell parameters for the X-ray experimental structure, the optimized Bulk Model structure and the optimized Isol Model one. The accuracy of the optimized Bulk Model internal degrees of freedom is remarkable as their deviation from the X-ray reference is about 1% or less.

Supplementary Figure 2.
First Coordination Shell Structural Features Scheme. Schematic representation of the F e 2+ (pink sphere) and N -ligands (blue sphere) geometry in [(tpa Ph )Fe] − . The symbols ri describes the distances between the central F e 2+ ion and the first coordination shell N-ligands. The α and β labels represent the NFeN angles. α is used for the angles involving only N atoms in the plane perpendicular to the pseudo 3-fold symmetry axis while β is used for angles involving the N atom along the pseudo 3-fold symmetry axis. The experimental structural features (Exp Model) are compared with the periodic DFT (Bulk Model) and non-periodic DFT (Isol Model) optimized values. The symbols ri describes the distances between the central F e 2+ ion and the first coordination shell N-ligands. The α and β labels represent the NFeN angles. α is used for the angles involving only N atoms in the plane perpendicular to the pseudo 3-fold symmetry axis while β is used for angles involving the N atom along the pseudo 3-fold symmetry axis.
The Bulk Model has been obtained through the optimization of both the ionic position inside the unit-cell and the lattice parameters. Supplementary Table 2 reports the comparison of the optimized lattice parameters with the experimental X-ray values. All the optimized lattice parameters show a very small deviation from the X-ray reference values, with errors of only a few percentage units. Moreover, it must be stressed that the X-ray structure has been collected at a temperature of 123 K [11], while our structures have been optimized at 0 K. Therefore, a shrinking of the simulated lattice parameters with respect to the X-ray reference is expected.
Normal modes and harmonic frequencies have been calculated for both the isolated and the periodic structure. The phonon density of states (DOSs) for both Isol Model and Bulk Model are displayed in Supplementary Fig. 3. Major differences between the two models are evident, especially in the low-energy part of the DOS. This clearly originates from the appearance of inter-molelcular vibrations. Interestingly, the Isol Model shows a non zero DOS at very low energy values. Further differences come from the presence in the periodic cell of the two Na + ions coordinated by three dimethoxyethanes each, which clearly contribute to the global DOS. Moreover, inter-molecular interactions are expected to introduce further modulation on the DOSs shape. CASSCF calculations have been performed in order to determine the [(tpa Ph )Fe] − zero-field splitting. The calculation of all the magnetic properties has been performed for a single isolated molecule as periodic boundary condition are not implemented for correlated calculations. For instance, for the Bulk Model only one of the two [(tpa Ph )Fe] − molecules in the primitive cell has been characterized. This does not represent a problem as they are symmetry-related inside the crystal. However, the elimination of the explicit surrounding of the SMM during the anisotropy calculation poses a serious limitation to the model as the crystal is made of charged molecules. The calculation of magnetic properties in SMMs is usually carried for isolated molecules and, therefore, electrostatic long-range effects originating from the presence of the crystal environment are usually neglected. However, the importance of second coordination shell effects on the magnetic structure of SMMs has been reported [12]. In order to check if the system under study shows these kind of features we introduced a point charges model to describe the charge density of the periodic crystal. Various point charge methods have been tested: Mulliken, density-derived point charges (DDAPC) [13], RESP [14] and CHELPG [15]. The first two methods rely on the projection of the DFT electronic structure on atomic localised functions, while the last two are based on a fitting of the electrostatic potential outside the van der Waals volume of the molecule coming from the DFT charge distribution. Although all these parametrisation schemes give similar results, final CASSCF calculations have been done with RESP point charges as this procedure is the one usually employed for force-field development and thus much more tested and robust. Supplementary Fig. 4 shows results about the Stark effect on the zero-field split (ZFS) originating from charged co-precipitate molecules. Results coming from the introduction in the CASSCF simulation cell of one of the nearest co-precipitate molecule (red bars) show the dramatic impact of the Stark effect on both the ZFS and its derivatives, the last one being calculated along a low-energy internal vibrational mode. The substitution of the explicit co-precipitate molecule with its fitted RESP point charges reproduces very nicely the results. It also must be stressed out that the explicit molecule model is not necessarily the most accurate one. Indeed the co-precipitate charged molecules is treated at the HF theory level, which is known to overestimate both the dipole moment and the polarisability [16]. Accordingly, DFT GGA fitted RESP point charges predict a smaller Stark effect. Calculations with the explicit co-precipitate molecule without the Na ions have also been carried out. The quenching of the net charge shows that dipole and polarisability effects have only a small impact on the magnetism of the SMM. In order to mimic the periodic surrounding, the single explicit SMM has been embedded in a 21×21×21 lattice of point charges. The model that includes the point charge correction is called Bulk PC Model. This model shares the same structure and phonon structure of the Bulk Model.

Supplementary
Before proceeding to the discussion of the D tensor calculations, a comment on the validity of the spin Hamiltonian itself is required. In a completely general fashion, the spin Hamiltonian formalism can be applied to describe the splitting of the spin energy levels as long as the wave-functions of the electronic states considered [see Supplementary Eq. (43)] possesses a well-defined S 2 expectation value, i.e. there is no spin contamination. Technically this becomes an issue only after the inclusion of the spin-orbit interaction, which may lead to a ground state multiplet with nonvanishing orbital angular momentum. In general Jahn-Teller activity removes the problem, as in the present case, where the global C 3 molecular symmetry is relaxed to a lower C 1 symmetry, as can also be seen from Supplementary  Table 1. However, the stiffness of the first coordination shell makes the Jahn-Teller distortions not fully effective in removing the degeneracy of the first two excited states and the final low-energy-lying spin-orbit-corrected wave-functions are indeed a superposition of the two different eigenstates of the Born-Oppenheimer Hamiltonian. However, in this case the two solutions have both an S=2 multiplicity and, therefore, this situation does not impose any formal restriction on the use of the spin Hamiltonian. It must also be stressed that the spin Hamiltonian formalism has also been already used before to study [(tpa Ph )Fe] − and other synthetic analogs [11,17]. Moreover, it is worth mentioning that our framework does not necessarily involve the use of the spin Hamiltonian. Indeed, one can directly use the expression (43) to evaluate the matrix elements in the first term of Supplementary Eq. (45), instead of passing by the fitting of the spin Hamiltonian. Supplementary Table 3 reports the calculated anisotropy parameters for all the models equilibrium structures. The D axial anisotropy constant calculated retaining the X-ray structure (X-ray model), without any optimisation, is found very close to the experimental one demonstrating the ability of CASSCF to handle this problem, as also previously noted by Atanasov et al. [17]. By looking at the various models it is possible to observe a small deviation from the X-ray value. Bulk Model is the one, which better reproduces the experimental value, according to its higher level of environment quality. Interestingly, the effect of the point charge corrections is almost negligible. Although preliminary calculations predict an important effect of an externally applied electrostatic field on the anisotropy, here the local field is everywhere zero due to the crystal symmetry and thus ineffective.
Although the calculated axial anisotropy parameters are in good agreement with the experimental values, the reproduction of the entire ground state S = 2 multiplet spectrum is less satisfactory. Indeed, energy levels of the second pseudo doublet, corresponding to DS 2 − D(S − 1) 2 = D(2 2 − 1 2 ) = −3D, are overestimated by about 15 cm −1 with respect to the CASSCF calculated ones, while the energy of the highest energy spin level,corresponding to DS 2 − D(S − 2) 2 = D(2 2 − 0 2 ) = −4D, is underestimated by about the same amount. The breaking down of the fit to a spin Hamiltonian containing only second order Stevens terms is in line with the high anisotropy of the molecule. Indeed, the large ZFS observed for this SMM originates from the presence of two almost degenerate electronic energy levels with the same spin multiplicity, which mix to create the S = 2 ground state spin multiplet once the SO interaction is turned on [17]. However, as the ground state multiplet does not suffer from S-mixing effects the quality of its description could be improved by employing higher-order terms in the spin Hamiltonian. Indeed, tests along this direction showed an increased quality of the CASSCF spin levels fitting when fourth order termŝ O kq ( S) are added to the spin Hamiltonian. However, the introduction of these terms does not change the results discussed in the next sections, except for an overall speeding up of the relaxation process thanks to new spin-phonon coupling terms. For this reason all the discussion on the spin dynamics is done over the data obtained with a spin HamiltonianĤ S = ij D ijŜiŜj without any loss of generality.
First and second order derivatives of the D tensor, with respect to the phonon variables, have been computed numerically by finite differences with a step (δq α ) of 0.2 normal mode units, according to the relations: where D ij and D ij (±δq α ) represents the ij element of the D tensor at the equilibrium geometry and at the distorted one, respectively. These two derivatives have also been extracted by interpolation of D ij elements as function of the displacement with a quadratic function f ij (q α ) = 1 2 ∂ 2 Dij ∂q 2 α q 2 α + ∂Dij ∂qα q α + k, where k is just an offset constant. The numerical value of both first and the second order derivative of the D ij elements, calculated with the two different methods, coincides for all the available digits for each phonon, assessing the quality of the calculated parameters. Supplementary Fig. 5 reports the D ij values as function of the displacement q α for the first normal mode (36 cm −1 ) of the Bulk Model and the corresponding fitted f ij (q α ) functions. The best fitted values obtained from the optimization of the f ij (q α ) functions, corresponding to the first and second order derivatives of the D ij parameters with respect to the q α normal mode with vibrational frequency 36 cm −1 , are reported in Supplementary Table 4 The determination of the second order derivatives of D with respect to a single normal modeq α also offers the possibility to calculate the average effect of the phonon bath dynamics on the spin Hamiltonian. This effect has been proposed by Pederson et al. [18] and very recently by Moreno et al. [19] to be at the origin of a spin Hamiltonian parameters renormalization effect. By assuming the spin-phonon coupling Hamiltonian to beĤ it is possible to evaluate the renormalization of the spin Hamiltonian parameters due to the average effect of the bath and defining a newĤ S asĤ Only the second order derivatives of the spin Hamiltonian contribute to the trace over the bath degrees of freedom appearing in Supplementary Eq. (48), leading tô According to Supplementary Table 4  The magnetisation dynamics has been studied by solving the Redfield equation for the diagonal elements of the spin reduced density matrix: where ρ S ab = a|ρ S |b , V α ab = a| ∂ĤS ∂qα |b , |a is the eigenfunction of the spin system described by the Hamiltonian H S with energy E a and ω ab = (E a − E b )/ . Thus in Supplementary Eq. (50) |V α ac | 2 G(ω ca , ω α ) represents the kinetic rate of population transfer between the eigenstates |a and |c .
From Supplementary Eq. (50) it is possible to calculate the evolution in time of the eigenstates population, which in turn can be used to compute the evolution in time of the magnetization expectation value as If we now choose of solving only the diagonal part of the Redfield Eqs. (50), we will restrict our study to the dynamics of pure states. In this case, the starting spin density matrixρ S (t = 0) must be built from one of the eigenvectors of the spin HamiltonianĤ S in order for the off-diagonal terms ofρ S to be zero at all times. The results displayed in the main text has been obtained by choosingρ S (t = 0) = |0 0|, where |0 corresponds to the lowest in energy eigenket ofĤ S . One can visualize this in silico experiment as the preparation of the quantum state by thermalizing the system at a temperature where only the ground state is populated and its afterward quench by abruptly increasing the temperature to the desired value. The evolution in time of the z-component of M is showed by the red line in Supplementary Fig. 6. The expectation value of M z at t = 0 is slightly lower than 2 due to the presence of the rhombic term inĤ S , which mixes theŜ z components in theĤ S eigenkets and lowers the magnetization expectation value. This effect is reduced by the presence of the external magnetic field (1500 Oe) present in the spin Hamiltonian but it can never be completely cancelled out. The time profile of M z is purely exponential and it has been used to extract the relaxation time, τ , from the fitting of M z (t) = (M z (0) − M z (∞))e −t/τ + M z (∞). In order to investigate the dependence of the spin relaxation time from the starting condition we also solved the full Redfield equation (52). This equation describes the evolution in time of all theρ S matrix elements and makes possible to study the evolution of arbitrary starting states. where By Considering the dynamics for only the diagonal part ofρ S one gets back the expression in Supplementary Eq. (50), i.e. R α aa,bb → 2M α a,b . The blue line of Supplementary Fig. 6 displays the dynamics of M z when the density matrix is initialised as pure maximally polarisedŜ z eigenstate. This scenario can be obtained in experiments by a preparation of the quantum state in an extremely high external magnetic field B followed by a quench due to the reduction of B to the desired value.
In this situation the exponential decay is superimposed to an oscillating function arising from the phase dynamics. This fast (∼ THz) wave function oscillations are the analog of a particle-through-barrier tunnelling effect. Indeed, if a generic quantum system is not initialised as an eigenstate, it will oscillate coherently between its eigenstates or part of them. This is exactly what happens to the spin system as long as it is not interacting with a bath. Clearly such oscillations are not part of a relaxation mechanism as the phase state dynamics is not damped. As soon as the interaction with a bath is switched on, every elements of the density matrix will start decaying. This effect must not be confused with the so called quantum tunnelling relaxation, which is instead a direct relaxation mechanism between two degenerate H S eigenstates due to inter molecule dipolar interactions. Quantum tunnelling relaxation is not covered by our simulations as only one SMM has been considered.
Even though the starting magnetisation is different in the two cases, the exponential-decay characteristic time is the same. This happens because the nature of the equilibrium driving force is unchanged in the two experiments and the population decaying part of the equation remains the same. For this reason all the calculations have been done by utilising only the diagonal elements of Supplementary Eq. (52), without any loss of generality. SUPPLEMENTARY NOTE 6. CP2K INPUT.
Here it follows the input file it has been used for periodic-DFT calculations with the software cp2k. As output of the vibrational analysis kind of calculation one gets frequencies and normal modes of the simulation cell in the molden format. Starting from these values is it possible to calculate the cartesian displacements associated to each normal mode and use them to perturb the Γ-point optimized geometry of a given molecule. This new set of distorted geometries has then be used as input for the ORCA software for the magnetic properties calculation.