Vibronic effects on the quantum tunnelling of magnetisation in Kramers single-molecule magnets

Single-molecule magnets are among the most promising platforms for achieving molecular-scale data storage and processing. Their magnetisation dynamics are determined by the interplay between electronic and vibrational degrees of freedom, which can couple coherently, leading to complex vibronic dynamics. Building on an ab initio description of the electronic and vibrational Hamiltonians, we formulate a non-perturbative vibronic model of the low-energy magnetic degrees of freedom in monometallic single-molecule magnets. Describing their low-temperature magnetism in terms of magnetic polarons, we are able to quantify the vibronic contribution to the quantum tunnelling of the magnetisation, a process that is commonly assumed to be independent of spin-phonon coupling. We find that the formation of magnetic polarons lowers the tunnelling probability in both amorphous and crystalline systems by stabilising the low-lying spin states. This work, thus, shows that spin-phonon coupling subtly influences magnetic relaxation in single-molecule magnets even at extremely low temperatures where no vibrational excitations are present.


I. INTRODUCTION
Single-molecule magnets (SMMs) hold the potential for realising high-density data storage and quantum information processing [1][2][3][4].These molecules exhibit a ground state comprising two states characterised by a large magnetic moment with opposite orientation, which represents an ideal platform for storing digital data.Slow reorientation of this magnetic moment results in magnetic hysteresis at the singlemolecule level at sufficiently low temperatures [5].The main obstacle to extending this behaviour to room temperature is the coupling of the magnetic degrees of freedom to molecular and lattice vibrations, often referred to as spin-phonon coupling [6].Thermal excitation of the molecular vibrations cause transitions between different magnetic states, ultimately leading to a complete loss of magnetisation.Advances in design, synthesis and characterisation of SMMs have shed light on the microscopic mechanisms underlying their desirable magnetic properties, and have allowed extending the nanomagnet behaviour to increasingly higher temperatures [7][8][9].
The mechanism responsible for magnetic relaxation in SMMs strongly depends on temperature.At higher temperatures, relaxation is driven by one (Orbach) and two (Raman) phonon transitions between magnetic sublevels [10].When temperatures approach absolute zero, all vibrations are predominantly found in their ground state.Thus, both Orbach and Raman transitions become negligible and the dominant mechanism is quantum tunnelling of the magnetisation (QTM) [11,12].This mechanism originates from a coherent coupling between the two magnetic ground states, which leads to the opening of a tunnelling gap.The tunnel coupling allows population to redistribute between states of opposite magnetisation, and thus facilitates magnetic reorientation.
While the role of vibrations in high-temperature magnetic relaxation is well understood in terms of weak-coupling rate equations for the electronic populations [13][14][15][16], the connection between QTM and spin-phonon coupling is still largely unexplored.Some analyses have looked at the influence of vibrations on QTM in integer-spin SMMs, where a model spin system was used to show that spin-phonon coupling could open a tunneling gap [17,18].However, QTM remains more elusive to grasp in half-integer spin systems, such as monometallic Dy(III) SMMs.In this case, a magnetic field is needed to break the time-reversal symmetry of the molecular Hamiltonian and lift the degeneracy of the ground doublet, as a consequence of Kramers theorem [19].This magnetic field can be provided by hyperfine interaction with nuclear spins or by dipolar coupling to other SMMs; both these effects have been shown to affect tunnelling behaviour [20][21][22][23][24][25][26][27].Once the tunnelling gap is opened by a magnetic field, molecular vibrations can in principle affect its magnitude in a nontrivial way (Fig. 1a).In a recent work, Ortu et al. analysed the magnetic hysteresis of a series of Dy(III) SMMs, suggesting that QTM efficiency correlates with molecular flexibility [23].In another work, hyperfine coupling was proposed to assists QTM by facilitating the interaction between molecular vibrations and spin sublevels [28].However, a clear and unambiguous demonstration of the influence of the spin-phonon coupling on QTM beyond toy-model approaches is still lacking to this date.A reason for this shortfall is found in the common wisdom that vibrations only cause transitions between electronic states when thermally excited, and therefore are unable to influence magnetic relaxation when thermal energy is much lower than their frequency.
In this work we present a theoretical analysis of the effect of molecular vibrations on the tunnelling dynamics in two prototypical Dy(III) SMMs, [Dy(Cp ttt ) 2 ] + [7] and [Dy(bbpen)Br] [29] (Fig. 1b).Our approach is based on a fully ab initio description of the SMM vibrational environment and accounts for the spin-phonon coupling in a non perturbative way.In this aspect, this work represents a step forward compared to previous theoretical analyses, which relied on a simplified description of phonons as small rotational displacements of the magnetic anisotropy axis and on a standard weak-coupling master equation approach [30].After deriving an effective lowenergy model for the relevant vibronic degrees of freedom based on a polaron approach [31], we demonstrate that vibrations can either enhance or reduce the quantum tunnelling gap, depending on the orientation of the magnetic field relative to the main anisotropy axis of the SMM.Lastly, we show that different vibrational modes can have competing effects on QTM; depending on how vibrations impact the axiality of the lowest energy magnetic doublet, they can lead to either a decrease or an increase of the tunnelling probability.While identifying vibrations that selectively tune QTM through chemical design of new SMMs goes beyond the scope of this work, our improved description of vibronic QTM provides a useful framework to articulate further studies in that direction.

A. Ab initio simulations
In this work we investigate two representative examples of Dy(III) SMMs and explore both amorphous and crystalline phonon environments.The first compound is [Dy(Cp ttt ) 2 ] + , shown in Fig. 1b, top [7].It consists of a dysprosium ion Dy(III) enclosed between two negatively charged cyclopentadienyl rings with tert-butyl groups at positions 1, 2 and 4 (Cp ttt ).The crystal field generated by the axial ligands makes the states with larger angular momentum be energetically favourable, resulting in the energy level diagram sketched in Fig. 1a.The energy barrier separating the two degenerate ground states results in magnetic hysteresis, which was observed up to T = 60 K [7].
To single out the contribution of molecular vibrations, we focus on a magnetically diluted sample in a frozen solution of dichloromethane (DCM).Thus, our computational model consists of a solvated [Dy(Cp ttt ) 2 ] + cation (Fig. 1b, top), which provides a realistic description of the low-frequency vibrational environment, comprised of pseudo-acoustic vibrational modes (Supplementary Note 1).These constitute the basis to consider further contributions of dipolar and hyperfine interactions to QTM.
Once the equilibrium geometry and vibrational modes of the solvated SMM (which are in general combinations of molecular and solvent vibrations) are obtained at the densityfunctional level of theory, we proceed to determine the equilibrium electronic structure via complete active space selfconsistent field spin-orbit (CASSCF-SO) calculations.The electronic structure is projected onto an effective crystal-field Hamiltonian.The spin-phonon couplings are obtained from a single CASSCF calculation by computing the analytic derivatives of the molecular Hamiltonian with respect to the nuclear coordinates [15].Further details can be found in the Methods section.
The second compound considered in this work is the highly stable [Dy(bbpen)Br] (H 2 bbpen= N, N ′ -bis(2hydroxybenzyl)-N, N ′ -bis(2-methylpyridyl)ethylenediamine), shown in Fig. 1b, bottom [29].It consists of a Dy(III) ion with pentagonal bipyramidal local geometry, with four N and one Br atom coordinating equatorially.Two axially coordinating O atoms give rise to strong easy-axis magnetic anisotropy.The effective barrier for magnetic reversal is around 1,000 K and magnetic hysteresis was observed up to 14 K [29].The small size of the unit cell and the relatively high-symmetry space group (C222 1 ) make this system amenable for spin-phonon coupling calculations in a crystalline environment.The primitive unit cell, consisting of two symmetry-related replicas of [Dy(bbpen)Br], was optimised at the density functional level of theory, and phonons were calculated using a 2 × 2 × 1 supercell expansion.The electronic structure of the Dy(III) centres was obtained with state-average CASSCF-SO and parametrised with a crystal field Hamiltonian.Spin-phonon couplings were obtained via the linear vibronic coupling model [15].A full account of these methods can be found in ref. [32].

B. Polaron model
The lowest-energy angular momentum multiplet of a Dy(III) SMM (J = 15/2) can be described by the ab initio vibronic Hamiltonian where E m denotes the energy of the m-th eigenstate |m⟩ of the crystal field Hamiltonian and Vj ⊗( b j + b † j ) represent the spinphonon coupling operators.The harmonic vibrational modes are described in terms of their bosonic annihilation (creation) operators b j ( b † j ) and frequencies ω j .In the absence of magnetic fields, the Hamiltonian (1) is symmetric under time reversal.This symmetry results in a two-fold degeneracy of the energy levels E m , whose corresponding eigenstates |m⟩ and | m⟩ form a time-reversal conjugate Kramers doublet.The degeneracy is lifted by introducing a magnetic field B, which couples to the electronic degrees of freedom via the Zeeman interaction ĤZee = µ B g J B • Ĵ, where g J is the Landé g-factor and Ĵ is the total angular momentum operator.To linear order in the magnetic field, each Kramers doublet splits into two energy levels E m ±∆ m /2 corresponding to the states where the energy splitting ∆ m and the mixing angles θ m and φ m are determined by the matrix elements of the Zeeman  6).Each spin state |1 ′ ± ⟩ is accompanied by a vibrational distortion (greatly exaggerated for visualisation), thus forming a magnetic polaron.Vibrational states |ν⟩ are now described in terms of harmonic displacements around the deformed structure, which depends on the state of the spin.Polarons provide an accurate physical picture when the spin-phonon coupling is strong and mostly modulates the energy of different spin states but not the coupling between them.
Hamiltonian on the subspace {|m⟩, | m⟩}.In addition to the intra-doublet mixing described by Eqs. ( 2) and (3), the Zeeman interaction also mixes Kramers doublets at different energies.The ground doublet acquires contributions from higherlying states These states no longer form a time-reversal conjugate doublet, meaning that the spin-phonon coupling can now contribute to transitions between them.
Since QTM is typically observed at much lower temperatures than the energy gap between the lowest and first excited doublets (which here is ≳ 600 K [7,29]) we focus on the perturbed ground doublet |1 ′ ± ⟩.Within this subspace, the Hamil-tonian Ĥ + ĤZee takes the form This Hamiltonian describes the interaction between vibrational modes and an effective spin one-half represented by the Pauli matrices describing the effect of the Zeeman interaction on the spin-phonon coupling.Due to the strong magnetic axiality of the SMM considered here, the longitudinal component of the spin-phonon coupling w z j dominates over the transverse part w x j , w y j .In this case, we can get a better physical picture of the system by transforming the Hamiltonian (5) to the polaron frame defined by the unitary operator which mixes electronic and vibrational degrees of freedom by displacing the mode operators by ξ ± j = (⟨1| Vj |1⟩ ∓ w z j )/ω j depending on the state of the effective spin one-half [31].The idea behind this transformation is to allow nuclei to relax around a new equilibrium geometry, which may be different for every spin state.This lowers the energy of the system and provides a good description of the vibronic eigenstates when the spin-phonon coupling is approximately diagonal in the spin basis (Fig. 1c).In the polaron frame, the longitudinal spin-phonon coupling is fully absorbed into the purely electronic part of the Hamiltonian, while the transverse components can be approximated by their thermal average over vibrations, neglecting their vanishingly small quantum fluctuations (Supplementary Note 2).After transforming back to the original frame, we are left with an effective spin one-half Hamiltonian with no residual spin-phonon coupling The set of Pauli matrices σ ′′ = Ŝ † (σ ′ ⊗ 1 vib ) Ŝ describe the two-level system formed by the magnetic polarons of the form Ŝ † |1 ′ ± ⟩|{ν j }⟩ vib , where {ν j } is a set of occupation numbers for the vibrational modes of the solvent-SMM system.These magnetic polarons can be thought as magnetic electronic states strongly coupled to a distortion of the molecular geometry.They inherit the magnetic properties of the corresponding electronic states, and can be seen as the molecular equivalent of the magnetic polarons observed in a range of magnetic materials [33][34][35].Polaron representations of vibronic systems have been employed in a wide variety of settings, ranging from spin-boson models [31,36] to photosynthetic complexes [37][38][39], to quantum dots [40][41][42], providing a convenient basis to describe the dynamics of quantum systems strongly coupled to a vibrational environment.These methods are particularly well suited for condensed matter systems where the electron-phonon coupling is strong but causes very slow transitions between different electronic states, allowing exact treatment of the pure-dephasing part of the electron-phonon coupling and renormalising the electronic parameters.For this reason, the polaron transformation is especially effective for describing our system (Supplementary Note 3).The most striking advantage of this approach is that the average effect of the spin-phonon coupling is included non-perturbatively into the electronic part of the Hamiltonian, leaving behind a vanishingly small residual spin-phonon coupling.
As a last step, we bring the Hamiltonian in Eq. ( 7) into a more familiar form by expressing it in terms of an effective gmatrix.We recall that the quantities ∆ 1 and w j depend linearly on the magnetic field B via the Zeeman Hamiltonian ĤZee .An additional dependence on the orientation of the magnetic field comes from the mixing angles θ 1 and φ 1 introduced in Eqs.
(2) and (3), appearing in the states |1 ± ⟩ used in the definition of w j .This further dependence is removed by transforming the Pauli operators back to the basis {|1⟩, | 1⟩} via a threedimensional rotation σ = R θ 1 ,φ 1 • σ ′′ .Finally, we obtain Ĥ(pol for appropriately defined electronic and single-mode vibronic g-matrices g el and g vib j .These are directly related to the electronic splitting term ∆ 1 and to the vibronic corrections described by w j in Eq. (7), respectively (see Supplementary Note 2 for a thorough derivation).The main advantage of representing the ground Kramers doublet with an effective spin one-half Hamiltonian is that it provides a conceptually simple foundation for studying low-temperature magnetic behaviour of the SMM, confining all microscopic details, including vibronic effects, to an effective g-matrix.

C. Vibronic modulation of the ground Kramers doublet
We begin by considering the influence of vibrations on the Zeeman splitting of the lowest Kramers doublet.The Zeeman splitting in absence of vibrations is simply given by In the presence of vibrations, the electronic g-matrix g el is modified by adding the vibronic correction ∑ j g vib j , resulting in the Zeeman splitting ∆ vib 1 .In Fig. 2a we show the Zeeman splittings as a function of the orientation of the magnetic field B for [Dy(Cp ttt ) 2 ] + , parametrised in terms of the polar angles (θ , φ ).Depending on the field orientation, vibrations can lead to either an increase or decrease of the Zeeman splitting.These changes seem rather small when compared to the largest electronic splitting, obtained when B is oriented along the z-axis (Fig. 1b), as expected for a system with easy-axis anisotropy.However, they become quite significant for field orientations close to the xy-plane, where the purely electronic splitting ∆ 1 becomes vanishingly small and ∆ vib 1 can be dominated by the vibronic contribution.This is clearly shown in Fig. 2b,c where we decompose the total field B = B int + B ext in a fixed internal component B int originating from dipolar and hyperfine interactions, responsible for opening a tunnelling gap, and an external part B ext which we sweep along a fixed direction across zero.When these fields lie in the plane perpendicular to the purely electronic easy axis, i.e. the hard plane, the vibronic splitting can be three orders of magnitude larger than the electronic one (Fig. 2b).The situation is reversed when the fields lie in the hard plane of the vibronic g-matrix (Fig. 2c).We note that this effect is specific to states with easy-axis magnetic anisotropy, however this is the defining feature of SMMs, such that our results should be generally applicable to all Kramers SMMs.In fact, we observe very similar results for [Dy(bbpen)Br] (Supplementary Note 4)., bottom) as a function of the orientation of the magnetic field, parametrised in terms of polar and azimuthal angles θ and φ .The polar angle θ is measured with respect to the axis joining the cyclopentadienyl centroids, corresponding approximately to the easy axis.The dashed (solid) line corresponds to the electronic (vibronic) hard plane.The magnitude of the magnetic field is fixed to 1 T. b, c, Electronic (dashed) and vibronic (solid) Zeeman splitting of the ground doublet as a function of the external field magnitude B ext in the presence of a transverse internal field B int = 1 mT calculated from Eq. ( 8).External and internal fields are perpendicular to each other and were both chosen to lie in the hard plane of either the electronic (b, purple) or vibronic (c, green) g-matrix.The orientation of the external (internal) field is shown for both cases as circles (crosses) in the inset in (a), with colors matching the ones in (b) and (c).

D. Internal fields and QTM probability
So far we have seen that spin-phonon coupling can either enhance or reduce the tunnelling gap in the presence of a magnetic field depending on its orientation.For this reason, it is not immediately clear whether its effects survive ensemble averaging in a collection of randomly oriented SMMs, such as for frozen solutions or polycrystalline samples considered in magnetometry experiments.In order to check this, let us consider an ideal field-dependent magnetisation mea-surement.When sweeping a magnetic field B ext at a constant rate from positive to negative values along a given direction, QTM is typically observed as a sharp step in the magnetisation of the sample when crossing the region around B ext = 0 [11,27].This sudden change of the magnetisation is due to a non-adiabatic spin-flip transition between the two lowest energy spin states, that occurs when traversing an avoided crossing (see diagram in Fig. 1a, right).The spin-flip probability is given by the celebrated Landau-Zener expression [43][44][45][46][47][48], which in our case takes the form where we have defined v = µ B dB ext /dt •g, and ∆ ⊥ is the component of ∆ = µ B B int • g perpendicular to v, while g denotes the total electronic-vibrational g-matrix appearing in Eq. ( 8) (see Supplementary Note 2 for a derivation of Eq. ( 9)).
In order to fully characterise the spin-flip process, we need to quantify the internal fields that cause QTM in Kramers SMMs, which originate from either dipolar or hyperfine interactions.In the following we focus on dipolar fields, since their effects can be observed at much higher temperatures than those required to witness hyperfine interactions (Supplementary Note 5).Samples studied in magnetometry experiments typically contain a macroscopic number of SMMs, each of which produces a microscopic dipole field.We estimate the combined effect of these microscopic dipoles in a [Dy(Cp ttt ) 2 ] + DCM frozen solution of by generating random spatial configurations of SMMs and calculating the resulting field at a specific point in space corresponding to a randomly selected SMM.We repeat this process 10,000 times to obtain the internal field distribution B int , as shown in Fig. 3a.The orientation of this field is random and its magnitude averages to 5.5 mT for a SMM concentration of 170 mM [7] (Supplementary Note 5).
In the case of the [Dy(bbpen)Br] molecular crystal, the effect of all Dy atoms within a 100 Å radius of a central magnetic centre was considered in a 5% Dy in Y diamagnetically diluted crystallite [29].Random Dy/Y subsitutions at different sites and random orientations of the magnetising field B ext were considered to mimic a powder sample, leading to the distribution shown in Fig. 3b with average magnitude 4.9 mT.
We then sample the distribution of internal fields to calculate the corresponding spin-flip probabilities for a randomly oriented SMM using Eq. ( 9).The effect of spin-phonon coupling on the spin-flip dynamics of an ensemble of SMMs is shown in Fig. 3c,d.The vibronic correction to the ground doublet g-matrix leads to a suppression of spin-flip events (orange) compared to a purely electronic model (blue).Despite the significant overlap between the two distributions, spinphonon coupling results in a ∼30% drop of average spin-flip probabilities, represented by the dashed lines in Fig. 3c,d.The vibronic suppression of QTM can be intuitively understood in terms of the polaron energy landscape sketched in Fig. 1c: strong coupling between spin degrees of freedom and molecular distortions can stabilise spin states, introducing a vibrational energy cost for spin reversal; i.e. flipping a spin requires reorganisation of the molecular structure.From Fig. 3c,d, we also note that crystalline [Dy(bbpen)Br] exhibits much larger QTM than [Dy(Cp ttt ) 2 ] + in frozen solution.This can be understood in terms of the different microscopic dipole fields in the two systems.In Supplementary Note 5 we show that B int is perfectly isotropic in a frozen solution.On the contrary, due to the symmetry of the [Dy(bbpen)Br] molecular crystal, the component of the internal field along the intra-unit cell Dy-Dy direction survives orientational averaging, resulting in an average transverse component of 1.2 mT (Supplementary Note 5).

III. DISCUSSION
As shown above, the combined effect of all vibrations in a randomly oriented ensemble of SMMs is to reduce QTM.However, not all vibrations contribute to the same extent.Based on the polaron model introduced above, vibrations with large spin-phonon coupling and low frequency have a larger impact on the magnetic properties of the ground Kramers doublet.This can be seen from Eq. ( 7), where the vibronic correction to the effective ground Kramers Hamiltonian is weighted by the factor ⟨1| Vj |1⟩/ω j .Another property of vibrations that can influence QTM is their symmetry.In monometallic SMMs, QTM has generally been correlated with a reduction of axial symmetry, either by the presence of flexible ligands or by transverse magnetic fields.Since we are interested in symmetry only as long as it influences magnetism, it is useful to introduce a measure of axiality on the g-matrix, such as where ∥ • ∥ denotes the Frobenius norm.This measure yields 1 for perfect easy-axis anisotropy, 1/2 for an easy-plane system, and 0 for the perfectly isotropic case.The axiality of an individual vibrational mode can be quantified as A j = A(g el +g vib j ) by building a single-mode vibronic g-matrix, analogous to the multi-mode one introduced in Eq. ( 8).We might be tempted to intuitively conclude that polaron formation always increases the axiality with respect to its electronic value A el = A(g el ), given that the collective effect of the spin-phonon coupling is to reduce QTM.However, when considered individually, some vibrations can have the opposite effect of effectively reducing the magnetic axiality.
In order to see how axiality correlates to QTM, we calculate the single-mode spin-flip probabilities ⟨P j ⟩.These are obtained by replacing the multi-mode vibronic g-matrix in Eq. ( 8) with the single-mode one g el + g vib j , and following the same procedure detailed in Supplementary Note 2. The single-mode contribution to the spin-flip probability unambiguously correlates with mode axiality, as shown in Fig. 4a for [Dy(Cp ttt ) 2 ] + ; the correlation is even starker for crystalline [Dy(bbpen)Br] (Fig. 4c).Vibrational modes that lead to a larger QTM probability are likely to reduce the magnetic axiality (top-left sector).Vice versa, those vibrational modes that enhance axiality also suppress QTM (bottom-right sector).
As a first step towards uncovering the microscopic basis of this unexpected behaviour, we single out the vibrational modes that have the largest impact on magnetic axiality in both directions.These vibrational modes, labelled A, B for [Dy(Cp ttt ) 2 ] + and C, D for [Dy(bbpen)Br], represent a range of qualitatively distinct vibrations, as can be observed in Fig. 4b,d.In the case of [Dy(Cp ttt ) 2 ] + , mode A is mainly localised on one of the Cp ttt ligands and features atomic displacements predominantly perpendicular to the easy axis.Mode B, on the other hand, involves axial distortions of the Cp rings and, to a lesser extent, rotations of the methyl groups.Thus, it makes sense intuitively that A would lead to an increased QTM probability, while the opposite is true for B, as observed in Fig. 4a.However, the connection between the magnetic axiality defined in Eq. ( 10) and vibrational motion is not always straightforward.In the case of [Dy(bbpen)Br], mode C mainly involves a tilt of the two equatorial pyridyl groups.This movement disrupts axiality and enhances QTM.On the other hand, mode D features equatorial motion of the first coordination sphere of the Dy(III) ion, involving movement of Br and Dy itself in the hard plane.However, this vibrational mode induces a suppression of QTM, as seen in Fig. 4c, rather than in increase, as would be expected based on the above symmetry arguments.This shows that ∆A j does not necessarily correlate to atomic motions, but can be a useful proxy for determining a given vibration's contribution to the QTM probability.In fact, the correlation between the two quantities can be rationalised with the help of the simple toy model presented in Supplementary Note 6.Nonetheless, we note that the out-of-phase motion of the equatorial pyridyl groups in D preserves axiality and could contribute to its efficiency at suppressing QTM.It is also worth noting that Briganti et al. recently demonstrated that motion of atoms beyond the first coordination sphere of the central Dy(III) ion can greatly influence spin dynamics in the Raman regime through bond polarisation effects [14].Performing a similar electrostatic analysis in the context of our polaron model is beyond the scope of this work; however, it represents an interesting direction for further investigations elucidating the role of vibrations on QTM.
In conclusion, we have presented a detailed description of the effect of molecular and solvent vibrations on the quantum tunnelling between low-energy spin states in two different single-ion Dy(III) SMMs, corresponding to amorphous and crystalline environments.Our theoretical results, based on an ab initio approach, are complemented by a polaron treatment of the relevant vibronic degrees of freedom, which does not suffer from any weak spin-phonon coupling assumption and is therefore well-suited to other strong coupling scenarios.We have been able to derive a non-perturbative vibronic correction to the effective g-matrix of the lowest-energy Kramers doublet, which we have used as a basis to determine the tunnelling dynamics in an idealised magnetic field sweep experiment, building on Landau-Zener theory.This has allowed us to formulate the observation that spin-phonon coupling does have an influence on QTM, albeit a subtle one (∼ 30%), as opposed to the widespread belief that magnetic tunnelling is not influenced by vibrations since it only becomes effective at low temperatures.This effect is rooted in the formation of magnetic polarons, which results in a redefinition of the magnetic anisotropy of the ground Kramers doublet.Our theoretical treatment is fully ab initio and represents a signifi-cant improvement over other theoretical descriptions of QTM which rely on weak coupling assumptions.Lastly, we observe that specific vibrational modes can either enhance or suppress QTM.This behaviour correlates to the magnetic axiality of each mode, which can be used as a proxy for determining whether a specific vibration enhances or hinders tunnelling.Our analysis suggests that there may be a positive side to spinphonon coupling in QTM.Enhancing the coupling to specific vibrations via appropriate chemical design while keeping detrimental vibrations under control, could in principle increase magnetic axiality and thus suppress QTM even further.However, translating this observation into clear-cut chemical design guidelines remains an open question, that requires the analysis of other molecular systems.As ab initio spin-phonon coupling calculations become more accessible, the approach presented here can be applied to the study of vibronic QTM in other SMMs, and thus represents a valuable tool for understanding the role of vibrations in low-temperature magnetic relaxation.

METHODS
The ab initio model of the DCM-solvated [Dy(Cp ttt ) 2 ] + molecule is constructed using a multi-layer approach.During geometry optimisation and frequency calculation the system is partitioned into two layers following the ONIOM scheme [49].The high-level layer, consisting of the SMM itself and the first solvation shell of 26 DCM molecules, is described by Density Functional Theory (DFT) while the outer bulk of the DCM ball constitutes the low-level layer modelled by the semi-empirical PM6 method.All DFT calculations are carried out using the pure PBE exchange-correlation functional [50] with Grimme's D3 dispersion correction.Dysprosium is replaced by its diamagnetic analogue yttrium for which the Stuttgart RSC 1997 ECP basis is employed [51].Cp ring carbons directly coordinated to the central ion are equipped with Dunning's correlation consistent triple-zeta polarised cc-pVTZ basis set and all remaining atoms with its double-zeta analogue cc-pVDZ [52].Subsequently, the electronic spin states and spin-phonon coupling parameters are calculated at the CASSCF-SO level explicitly accounting for the strong static correlation present in the f-shell of Dy(III) ions.At this level, environmental effects are treated using an electrostatic point charge representation of all DCM atoms.All DFT/PM6 calculations are carried out with GAUSSIAN version 9 revision D.01 [53] and the CASSCF calculations are carried out with OPENMOLCAS version 21.06 [54].
The starting [Dy(Cp ttt ) 2 ] + solvated system was obtained using the solvate program belonging to the AmberTool suite of packages, with box as method and CHCL3BOX as solvent model.Chloroform molecules were subsequently converted to DCM.From this large system, only molecules falling within 9 Å from the central metal atom are considered from now on.The initial disordered system of 160 DCM molecules packed around the [Dy(Cp ttt ) 2 ] + crystal structure [7] is pre-optimised in steps, starting by only optimising the high-level layer atoms and freezing the rest of the system.The low-layer atoms are pre-optimised along the same lines starting with DCM molecules closest to the SMM and working in shells towards the outside.Subsequently, the whole system is geometry optimised until RMS (maximum) values in force and displacement corresponding to 0.000 45 au (0.0003 au) and 0.0018 au (0.0012 au) are reached, respectively.After adjusting the isotopic mass of yttrium to that of dysprosium m Dy = 162.5 u, vibrational normal modes and frequencies of the entire molecular aggregate are computed within the harmonic approximation.
Electrostatic atomic point charge representations of the environment DCM molecules are evaluated for each isolated solvent molecule independently at the DFT level of theory employing the CHarges from ELectrostatic Potentials using a Grid-based (ChelpG) method [55], which serve as a classical model of environmental effects in the subsequent CASSCF calculations.
The evaluation of equilibrium electronic states and spinphonon coupling parameters is carried out at the CASSCF level including scalar relativistic effects using the secondorder Douglas-Kroll Hamiltonian and spin-orbit coupling through the atomic mean field approximation implemented in the restricted active space state interaction approach [56,57].The dysprosium atom is equipped with the ANO-RCC-VTZP, the Cp ring carbons with the ANO-RCC-VDZP and the remaining atoms with the ANO-RCC-VDZ basis set [58].The resolution of the identity approximation with an on-the-fly acCD auxiliary basis is employed to handle the two-electron integrals [59].The active space of 9 electrons in 7 orbitals, spanned by 4f atomic orbitals, is employed in a state-average CASSCF calculation including the 18 lowest lying sextet roots which span the 6 H and 6 F atomic terms.
We use our own implementation of spin Hamiltonian parameter projection to obtain the crystal field parameters B q k entering the Hamiltonian describing the 6 H 15/2 ground state multiplet.Operator equivalent factors and Stevens operators are denoted by θ k and O q k ( Ĵ), where Ĵ = ( Ĵx , Ĵy , Ĵz ) are the angular momentum components.Spin-phonon coupling arises from changes to the Hamiltonian (11) due to slight distortions of the molecular geometry, parametrised as where X j denotes the dimensionless j-th normal coordinate of the molecular aggregate.The derivatives ∂ B q k /∂ X j are calculated using the Linear Vibronic Coupling (LVC) approach described in Ref. [15] based on the state-average CASSCF density-fitting gradients and non-adiabatic coupling involving all 18 sextet roots.Finally, we express the dimensionless normal coordinates in terms of bosonic creation and annihilation operators as Xj = ( b j + b † j )/ √ 2, which defines the system part of the spin-phonon coupling operators in Eq. ( 1) as DATA AVAILABILITY The data generated in this study have been deposited in the Figshare database and can be accessed at http://doi.org/10.48420/21892887[60].Source data for all figures are provided with this paper.

B. Polaron Hamiltonian for the ground doublet
Now that we have an approximate expression for the relevant electronic states, we reintroduce the spin-phonon coupling into the picture.First, we project the vibronic Hamiltonian (S1) onto the subspace spanned by |1 ′ ± ⟩, yielding On this basis, the purely electronic part ĤCF + ĤZee is diagonal with eigenvalues E 1 ± ∆ 1 /2, and the purely vibrational part is trivially unaffected.On the other hand, the spin-phonon couplings can be calculated to lowest order in the magnetic field strength B as where we have defined Ŵj = Vj Q1 ĤZee + ĤZee Q1 Vj (S13) and used the time-reversal invariance of the spin-phonon coupling operators to obtain The two states |1 ± ⟩ form a conjugate pair under time reversal, meaning that Θ|1 ± ⟩ = ∓e iα |1 ∓ ⟩ for some α ∈ R. Using the fact that for any two states ψ, ϕ, and for any operator Ô we have ⟨ψ| Ô|ϕ⟩ = ⟨ Θϕ| Θ Ô † Θ−1 | Θψ⟩, and recalling that the angular momentum operator is odd under time reversal, i.e.ΘĴ Θ−1 = − Ĵ, we can show that Keeping in mind these observations, and defining the vector we can rewrite the spin-phonon coupling operators in Eq. (S10) as where σ ′ is a vector whose entries are the Pauli matrices in the basis Plugging this back into Eq.(S10) and explicitly singling out the diagonal components of Ĥeff in the basis |1 ′ ± ⟩, we obtain At this point, we apply a unitary polaron transformation to the Hamiltonian (S16) where ξ s j = ⟨1| Vj |1⟩ − sw z j /ω j and is the bosonic displacement operator acting on mode j, i.e.D j (ξ ) b j D † j (ξ ) = b j − ξ .The Hamiltonian thus becomes The polaron transformation reabsorbes the diagonal component of the spin-phonon coupling (S15) proportional to w z j into the energy shifts ω j |ξ ± j | 2 , leaving a residual off-diagonal spin-phonon coupling proportional to w x j and w y j .Note that the polaron transformation exactly diagonalises the Hamiltonian (S10) if w x j = w y j = 0.In Supplementary Note 3, we argue in detail that in our case |w x j |, |w y j | ≪ |w z j | to a very good approximation.Based on this argument, we could decide to neglect the residual spin-phonon coupling in the polaron frame.The energies of the states belonging to the lowest doublet are shifted by a vibronic correction leading to a redefinition of the energy gap Although the off-diagonal components of the spin-phonon coupling w x j and w y j are several orders of magnitude smaller than the diagonal one w z j (see Supplementary Note 3), the sheer number of vibrational modes could still lead to an observable effect on the electronic degrees of freedom.We can estimate this effect by averaging the residual spin-phonon coupling over a thermal phonon distribution in the polaron frame.Making use of Eq. (S17), the off-diagonal coupling in Eq. (S19) can be written as Assuming the vibrations to be in a thermal state at temperature T in the polaron frame obtaining the average of Eq. (S23) reduces to calculating the dimensionless quantity which appears as a multiplicative rescaling factor for the off-diagonal couplings ⟨1 ∓ | Ŵj |1 ± ⟩.Note that, when neglecting second and higher order terms in the magnetic field, κ j does not show any dependence on temperature or on the magnetic field orientation via θ 1 and φ 1 .
After thermal averaging, the effective electronic Hamiltonian for the lowest energy doublet becomes where the energy of the lowest doublet is shifted by due to the spin-phonon coupling and to the thermal phonon energy.Eq. (S26) thus represents a refined description of the lowest effective spin-1/2 doublet in the presence of spin-phonon coupling.We can finally recast the Hamiltonian (S26) in terms of a g-matrix for an effective spin 1/2, similarly to what we did earlier in the case of no spin-phonon coupling.In order to do so, we first recall from Eq. (S6) and (S14) that the quantities ∆ 1 and (w x j , w y j , w z j ) appearing in Eq. (S26) depend on the magnetic field orientation via the states |1 ± ⟩, and on both orientation and intensity via ĤZee .We can get rid of the first dependence by expressing the Zeeman eigenstates |1 ± ⟩ in terms of the original crystal field eigenstates |1⟩, | 1⟩.For the spin-phonon coupling vector w j , we obtain (S28) where R(θ 1 , φ 1 ) is a rotation matrix.Similarly, the elctronic contribution ∆ 1 transforms as The Pauli spin operators need to be changed accordingly to σ = R(θ 1 , φ 1 ) T • σ ′ .Lastly, we single out explicitly the magnetic field dependence of Ŵj , defined in Eq. (S13), by introducing a three-component operator K j = ( Kx j , Ky j , Kz j ), such that Thus, the effective electronic Hamiltonian in Eq. (S26) can be finally rewritten as el + g vib • σ/2 (S31) where g (1) el is the electronic g-matrix defined in Eq. (S3), and is a vibronic correction.Note that this correction is non-perturbative in the spin-phonon coupling, despite only containing quadratic terms in Vj (recall that K j depends linearly on Vj ).The only approximations leading to Eq. (S31) are a linear perturbative expansion in the magnetic field B and neglecting quantum fluctuations of the off-diagonal spin-phonon coupling in the polaron frame, which is accounted for only via its thermal expectation value.This approximation relies on the fact that the off-diagonal couplings are much smaller than the diagonal spin-phonon coupling that is treated exactly by the polaron transformation (see Supplementary Note 3).

C. Landau-Zener probability
Let us consider a situation in which the magnetic field comprises a time-independent contribution arising from internal dipolar or hyperfine fields B int and a time dependent external field B ext (t).Let us fix the orientation of the external field and vary its magnitude at a constant rate, such that the field switches direction at t = 0.Under these circumstances, the Hamiltonian of Eq. (S31) becomes

FIG. 1 .
FIG. 1.Quantum tunnelling in Dy(III) single-molecule magnets.a, Typical energy level diagram of the lowest-energy J multiplet with angular momentum J = 15/2 in a Dy(III) single-molecule magnet (SMM), with degenerate doublets at energies E 1 , E 2 , etc. States are organised according to the expectation value of the total angular momentum along the magnetic anisotropy axis ⟨ Ĵz ⟩.Dipolar and hyperfine magnetic fields (B int ) can lift the degeneracy of the ground doublet and cause quantum tunnelling of the magnetisation (QTM), which results in avoided crossings when sweeping an external magnetic field B ext .Molecular vibrations can influence the magnitude of the energy splitting ∆ 1 .b, Top: Molecular structure of [Dy(Cp ttt ) 2 ] + surrounded by a dichloromethane (DCM) bath.Bottom: Structure of a [Dy(bbpen)Br] molecular crystal.Only the two SMMs in the primitive unit cell are shown; violet spheres represent Dy atoms at other lattice positions.Atoms are colour coded as follow: Dy (violet), Br (brown), Cl (green), O (red), N (cyan), C (grey), H (white).In both cases, z indicates the direction of the easy axis.c, Idea behind the polaron transformation Ŝ of Eq. (6).Each spin state |1 ′ ± ⟩ is accompanied by a vibrational distortion (greatly exaggerated for visualisation), thus forming a magnetic polaron.Vibrational states |ν⟩ are now described in terms of harmonic displacements around the deformed structure, which depends on the state of the spin.Polarons provide an accurate physical picture when the spin-phonon coupling is strong and mostly modulates the energy of different spin states but not the coupling between them.

FIG. 2 .
FIG. 2. Zeeman splitting of the ground Kramers doublet in[Dy(Cp ttt ) 2 ] + .a, Electronic ground doublet splitting (∆ 1 , top) and vibronic correction (∆ vib 1 − ∆ 1 , bottom) as a function of the orientation of the magnetic field, parametrised in terms of polar and azimuthal angles θ and φ .The polar angle θ is measured with respect to the axis joining the cyclopentadienyl centroids, corresponding approximately to the easy axis.The dashed (solid) line corresponds to the electronic (vibronic) hard plane.The magnitude of the magnetic field is fixed to 1 T. b, c, Electronic (dashed) and vibronic (solid) Zeeman splitting of the ground doublet as a function of the external field magnitude B ext in the presence of a transverse internal field B int = 1 mT calculated from Eq. (8).External and internal fields are perpendicular to each other and were both chosen to lie in the hard plane of either the electronic (b, purple) or vibronic (c, green) g-matrix.The orientation of the external (internal) field is shown for both cases as circles (crosses) in the inset in (a), with colors matching the ones in (b) and (c).

FIG. 3 .
FIG. 3. Internal fields and spin-flip probability.a, b, Distribution of internal field magnitudes B int experienced by a Dy centre due to the dipolar fields produced by surrounding Dy centres, magnetised by a randomly oriented external field B ext .For [Dy(Cp ttt ) 2 ] + (a), a uniform spatial distribution of 1,000 randomly oriented singlemolecule magnets (SMMs) around a central Dy(III) was assumed, corresponding to a 170 mM solution in dichloromethane.For the [Dy(bbpen)Br] molecular crystal (b), we considered the total dipolar field arising from all Dy centres within a 100 Å radius from a central Dy assuming 5% diamagnetic dilution.c, d, Distribution of electronic (blue) and vibronic (orange) Landau-Zener spin-flip probabilities P LZ , calculated for a randomly oriented SMM subjected to the dipolar fields shown above, assuming an external field sweep rate of 10 Oe/s.Average values are shown as dashed lines: (c) 0.0104 (blue) and 0.0074 (orange); (d) 0.903 (blue) and 0.618 (orange).All histograms are obtained from an ensemble of 10,000 random external field orientations and dipole arrangements.

FIG. 4 .
FIG. 4. Single-mode contributions to tunnelling of the magnetisation.a, c, Single-mode vibronic spin-flip probabilities plotted for each vibrational mode, shown as a function of the mode axiality ∆A j = A j − A el relative to the electronic axiality A el .The magnitude of the internal field is fixed to B int = 1 mT and the external field sweep rate is 10 Oe/s.The probabilities ⟨P j ⟩ are obtained by averaging over random orientations of external and internal fields.The color coding represents the spin-phonon coupling strength ∥ Vj ∥.Grey dashed lines corresponds to a purely electronic model.(a) and (c) correspond to amorphous [Dy(Cp ttt ) 2 ] + and crystalline [Dy(bbpen)Br].b, d Visual representation of the displacements induced by the vibrational modes indicated by arrows in (a) and (c) denoted by A, B, C, D; the corresponding vibrational frequencies are denoted by ω A ,ω B ,ω C ,ω D .