The mechanism of spin-phonon relaxation in endohedral metallofullerene single molecule magnets

Abstract

This study presents the investigation of spin-phonon coupling mechanisms in fullerene-based single-molecule magnets (SMMs) using ab initio CASSCF combined with DFT calculations. While lanthanide-based SMMs, particularly those with Dy^III ions, are known for their impressive blocking temperatures and relaxation barriers, endohedral metallofullerene (EMFs) offer a unique platform for housing low-coordinated lanthanides within rigid carbon cages. We have explored the spin dynamics in DyScS@C₈₂ exhibiting among the highest blocking temperature (T_B) reported. Through our computational analysis, we reveal that while the fullerene cage enhances crystal field splitting and provides structural stability without significantly contributing to spin-relaxation driving low-energy phonons, the internal ionic motion emerges as the primary factor controlling spin relaxation and limiting blocking temperature. This computational investigation into the spin dynamics of EMF-based SMMs provides key insights into their magnetic behaviour and suggests potential strategies for improving their performance towards futuristic SMMs.

---

Introduction

Due to their very high blocking temperatures and large thermal barriers to relaxation, Dy^III-based single-molecule magnets (SMMs) are of great interest.^1–6 Among many challenges in taking these molecules to end-user applications, enhancing the barrier height for magnetization reversal (U_eff) and blocking temperature below which the magnetisation is fully frozen (T_B) are considered crucial. Coupling magnetic lanthanide ions to strong axial ligand fields is one of the most promising approaches to producing high-performing SMMs.^1,3 This has produced SMMs with a barrier height of magnetisation reversal as large as 1500 cm⁻¹, achieving one of the aforementioned goals.^1,2 Although such large barriers have been achieved so far the T_B values remain modest at 80 K.¹ While theoretical advances have propelled the understanding of U_eff through various intuitive ligand design principles that have shaped the field of lanthanide SMMs over the past few decades, the comprehension of T_B, the most critical factor governing magnetic performance, remains largely unclear. This is evident from the fact that many molecules with high U_eff display low T_B, and vice versa, underscoring the need for a molecular-level understanding beyond simple electron density considerations (such as prolate versus oblate shapes). To achieve the next breakthrough in enhancing T_B values, it is essential to focus on this deeper molecular understanding. The T_B values and the corresponding relaxation times are strongly linked to the spin-phonon relaxation mechanism, an area of significant recent research aimed at providing insights to improve T_B.^1,2,7–9

The mechanism of spin-phonon relaxation has been explored so far in a handful of examples offering some guidelines to enhance the T_B values.^8,10,11 These guidelines include: (i) achieving strong crystal-field splitting of the Kramers doublets,^12,13 (ii) ensuring that the crystal field splitting exhibits minimal transverse anisotropy, which drives quantum tunnelling of magnetization (QTM),^14,15 (iii) reducing the vibrational density of states at resonant frequencies,¹⁴ (iv) minimizing low-energy vibrations,¹⁶ (v) utilizing rigid ligands that can isolate intramolecular motions from low-energy acoustic vibrations,^7,17 and (vi) employing ligands with donor atoms whose local charge remains stable despite local vibrations.⁸ Controlling these factors is very challenging, as crystal field parameters and reducing transverse anisotropy are correlated to the geometry and local point group of the molecules. Strong crystal field splitting and less transverse anisotropy can be achieved with low-coordinated lanthanides or with higher oxidation state,¹⁸ but synthesizing such molecules is challenging, and even if they are synthesised with very bulky ligands, they often encounter other weak non-covalent interactions such as Ln⋯H–C agostic interactions that facilitate transverse anisotropy as demonstrated in several cases.^2,19,20 Moreover, traditional organometallic complexes often have ligands with loosely bound atoms, like hydrogen, whose vibrations also contribute to spin-phonon relaxation.² Synthesis of the rigid lanthanide low-coordinate molecule without atoms such as –H are extremely challenging and has not yet been achieved.

However, in the endohedral metallofullerenes (EMFs) class of lanthanide-based SMMs, these conditions can be met easily.^21–26 The EMF-based SMMs feature rigid cages that can stabilize low-coordinated lanthanides inside their cage structure.^26,27 These compounds contain atoms or small clusters of atoms encapsulated by fullerenes, where cage shields the magnetic ion from the decoherence caused by external noise and hence stabilize atomic configurations that are not achievable in conventional molecules. These unique circumstances can produce very effective, precisely controlled SMMs. Although not all conditions are met, this type of molecule often satisfies conditions (i), (ii), (iv) (v) and (vi). This is demonstrated by the fact that lanthanide EMFs are reported as high-performing SMMs far more frequently than any other class of molecules (see Table S1†). Despite these key advantages, the T_B that has been reached with this class of molecule is relatively small compared to the best-in-class dysprosocenium class of molecules. This is linked to the spin-phonon relaxation mechanism, but a clear understanding of this process for this class of SMMs is still missing. Investigating the spin-phonon relaxation mechanism in these fullerene molecules could provide valuable insights, potentially leading to methods for suppressing or minimizing this effect.

Endohedral metallofullerenes can be broadly classified into two main categories. The first category consists of traditional metallofullerenes, where only metal atoms are encapsulated within the fullerene structure. In contrast, the second category, known as cluster fullerenes, encloses not only metal atoms but also non-metals such as nitrides, oxides, sulphides, carbides, and related cyanides within the fullerene shell.^22,28–37 Although among all fullerene based SMMs, Dy@C_s(6)-C₈₁N possess the largest T_B value,²⁵ among cluster fullerenes with single magnetic centre, DyScS@C₈₂ exhibit the one of the best T_B value of 7.3 K. Almost all reported EMFs exhibiting SMM characteristics within this range, suggesting a common spin-phonon relaxation mechanism in this class of molecules that may limit the achievement of higher T_B values. In this manuscript, we investigated the spin-phonon relaxation processes in one such cluster fullerene DyScS@C₈₂ using state-of-the-art ab initio spin relaxation simulations. To the best of our knowledge, our study marks the first implementation of ab initio spin-phonon relaxation calculations in the context of fullerene-based SMMs. Our goal is to (i) elucidate the mechanism of spin-phonon relaxation, (ii) compute relaxation times for comparison with experimental data, (iii) analyse the vibrational modes responsible for spin-phonon relaxation, and (iv) provide general guidance on how to achieve larger T_B values in this class of molecules.

In the original paper by Echegoyen and co-workers,³⁷ the synthesis and characterisation of two new dysprosium-containing mixed-metallic sulphide cluster fullerenes, DyScS@C_s(6)–C₈₂ and DyScS@C_3v(8)−C₈₂ were reported. These compounds were isolated and analysed using various techniques, including mass spectrometry, Vis-NIR spectroscopy, cyclic voltammetry, and single-crystal X-ray diffraction. For this manuscript, we chose DyScS@Cs(6)−C₈₂ for our calculations due to its higher U_eff barrier and longer relaxation times in the absence of external magnetic fields compared to DyScS@C_3v(8)−C₈₂. In the remainder of the manuscript, we will refer to DyScS@C_s(6)-C₈₂ simply as DyScS@C₈₂.

Computational details

To execute geometry optimizations, as well as simulations of Γ-point phonons, the CP2k program package^38–40 was used. The CP2k software uses a hybrid basis set technique called the Gaussian and Plane Wave Method (GPW),^39,40 in which the electronic charge density is represented by an auxiliary plane-wave basis set and Kohn–Sham orbitals are enlarged using contracted Gaussian-type orbitals (GTOs). A large cell of 20 ų containing a single fullerene unit was used for geometry optimization in absence of periodic boundary conditions. Utilising their relativistic norm-conserving pseudopotentials (Goedecker, Teter, and Hutter)⁴¹ optimised for the PBE functional.⁴² The Dy^III ion was replaced by the Y^III ion to avoid convergence issues with the DFT framework, however the atomic mass of Dy was used in order to create the vibrational modes. The DZVP-MOLOPT-GTH basis set (valence double-zeta (ζ) plus polarisation, molecularly optimised, Goedecker–Teter–Hutter) for all atoms (H, C) and the DZVP-MOLOPT-SR-GTH basis set (valence double-zeta (ζ) plus polarisation, molecularly optimised, short range Goedecker–Teter–Hutter) for Y, as implemented in the CP2k, were used for all calculations.^40,43 Furthermore, 1000 Ry was used as the energy cut-off for the plane wave basis set. For geometry optimization, a very tight force convergence criterion of 10⁻⁸ au and an SCF convergence criterion of 10⁻⁸ au for energy were employed (see input file in ESI†). Here it is important to discuss that the EMFs are known to exhibit a significant degree of structural disorder due to the dynamic nature of the encapsulated metal cluster and the flexibility of the fullerene cage. In our study, we have initiated the modeling from an optimized geometry to provide a well-defined and reproducible starting point for theoretical analysis. While this approach offers clarity and computational tractability, it does not fully capture the inherent disorder present in experimental systems. To better address these complexities, molecular dynamics (MD) simulations can be employed to explore the range of accessible conformations and mimic dynamic disorder.⁵ Such techniques although beyond the scope of this paper, offer a promising route to understand how structural fluctuations influence the magnetic properties of EMF-based single-molecule magnets.

The crystal field parameters for the geometries were computed using the ORCA 5.0 suite of software.⁴⁴ For the Dy^III ions, magnetic properties were derived from CASSCF calculations with an active space of seven 4f orbitals containing nine electrons (9,7), considering all solutions with multiplicities of 6. The basis sets DKH-SVP for C atoms, and SARC-DKH-TZVP for Dy atoms,⁴⁵ DKH-TZVP for Sc and S atoms were utilized.^46,47 Also, we have performed CASSCF/RASSI-SO⁴⁸/SINGLE_ANISO calculations using the Molcas 8.2 package in order to compute the static relaxation mechanism in DyScS@C₈₂ molecule and the {Dy–S–Sc}⁴⁺ fragment.^49,50 We have employed the ANO-RCC-TZVP basis set⁵¹ for Dy^III ions and ANO-RCC-DZVP for rest of the atoms. We have used 21 sextets in CASSCF and RASSI-SO calculations as demonstrated earlier in this class of compounds.^5,52,53 We have also used the Gaussian 16 package⁵⁴ for the Atoms in Molecules (AIM) analysis, employing the UB3LYP functional^55–57 with the SDD basis set⁵⁸ for Y, CSDZ^59,60 for Dy^III, and 6-31G* for the remaining atoms.⁶¹ Given the heavy computational requirements of spin-phonon calculations, due to the need to sample many geometries and electronic states, dynamical correlation was excluded from our study. However, previous investigations^18,52,62 support that this approximation maintains qualitative consistency in magnetic trends.

Calculations of the crystal field and spin-phonon coupling coefficients

The effective crystal-field Hamiltonian employed here correspond to

The parameters B^l_mare the CF Hamiltonian coefficients, and the operators Ô^l_m are tesseral functions of the total angular momentum operator J, with rank l and order m.⁶³ The lowest 2J+1 ab initio wave functions are matched one-to-one with the magnetic ground state of the ion to determine the CF Hamiltonian coefficients from first principles. To achieve this, the spin Hamiltonian |J̃m_J〉 is obtained by diagonalizing the operator Ĵ_z within the basis of the lowest 2J+1 ab initio wave functions.⁶⁴ Subsequently, the ab initio energy matrix is represented in this basis, and the parameters in eqn (1) are adjusted to reproduce the ab initio energy matrix elements.⁶⁴ These values enable the calculation of the magnetic anisotropy tensor across all temperatures and the magnetization in every spatial direction. The crystal field Hamiltonian, as described in eqn (1), represents the energy levels of the ground-state electronic multiplet for the equilibrium geometry of the molecule. However, due to thermal energy, the molecular geometry is constantly fluctuating. These fluctuations lead to the modulation of spin properties, namely spin-phonon coupling, causing transitions between different electronic states until a thermal equilibrium is achieved between the electronic states and the lattice. Therefore, it is necessary to describe the spin, the phonons, and their interaction in detail in order to provide a quantum mechanical explanation of spin-dynamics. To derive the spin-phonon coupling coefficients (∂Ĥ_s/∂Q_α), B^l_m that occur in eqn (1) are numerically differentiated with respect to the atomic displacements described by Q_α.

where N is the number of atoms in the unit cell, Q_α is the displacement vector connected to the α-phonon, and L_iα and ω_α are the eigenvectors of the Hessian matrix and the angular frequency of the phonon, respectively. By using numerical differentiation, the spin Hamiltonian's first-order derivatives (∂ĥ_S/∂X_i) with respect to the Cartesian degree of freedom X_i are determined. Four samples are taken between ±0.1 Å for each molecular degree of freedom.

Calculation of spin-relaxation

The total Hamiltonian operator is defined as:

Ĥ = Ĥ_S + Ĥ_Ph + Ĥ_S−Ph (3)

where the spin and phonon Hamiltonians are represented by the first two terms, respectively, and the coupling between these subsystems is represented by the third term. All-electron system's low-lying electronic states are described by the spin Hamiltonian. A straightforward sum of harmonic oscillators is used to simulate the phonon Hamiltonian. It is worth noting that while incorporating anharmonic interactions would provide a more complete picture, doing so requires computing linewidths from first-principles methods and perturbation theory— an approach that remains computationally expensive and not necessarily more accurate.^8,65–68

After obtaining the eigenstates, |a〉, and eigenvalues, E_a, of these operators, spin dynamics can be simulated by calculating the W_ab, or transition rate between distinct spin states.⁶⁹ In molecular Kramers systems with high magnetic anisotropy, one- and two-phonon processes contribute to spin relaxation. When one-phonon processes are taken into account, the transition rate between spin states, Ŵ^1−Ph_ba is

where the phrase (∂Ĥ₀/∂Q_α) indicates the strength of the coupling between spin and the α-phonon Q_α, and ℏω_ba = E_b – E_a. The G^1−Ph function reads

G^1−Ph(ω,ω_α) = δ(ω − ω_α)n̄_α + δ(ω + ω_α)(n̄_α+ 1) (5)

where the Bose–Einstein distribution accounting for the thermal population of phonons is represented by , the Boltzmann constant is represented by k_B, and the Dirac δ functions enforce energy conservation during the absorption and emission of phonon by the spin system, respectively, here δ functions has been approximated using Gaussian smearing. The finite differentiation is used to calculate lattice harmonic frequencies (ω_α/2π) and normal modes (Q_α) following geometry optimisation using DFT. The Orbach relaxation mechanism is explained by eqn (2), where the spin moves from the fully polarised state M_J = J to an excited state with an intermediate value of M_J before it can emit phonons again and return to M_J = −J. This process also occurs for states represented by the total angular momentum, J. A further route of relaxation to equilibrium is offered by two-phonon processes, leading to the Raman mechanism. Two-phonon spin-phonon transitions, or Ŵ^2−Ph_baare modelled aswhere T^αβ,±_ba isT^αβ,±_ba entails the simultaneous contribution of all spin states |c〉, often known as a virtual state. All two-phonon processes, such as the absorption of two phonons, the emission of two phonons, or the absorption of one phonon and the emission of a second, are taken into account by G^2−Ph. The latter method, which in this instance yields a G^2−Ph reading, determines the Raman relaxation rate.

G^2−Ph (ω,ω_α,ω_β) = δ(ω − ω_α + ω_β) n̄_α(n̄_β + 1) (8)

After computing each matrix element Ŵ^n−Ph_ba, τ⁻¹ can be anticipated by diagonalizing Ŵ^n−Ph_ba and determining the smallest nonzero eigenvalue. The Orbach contribution to the relaxation rate, τ⁻¹_Orbach, is obtained from the analysis of Ŵ^1−Ph, while the Raman contribution, τ⁻¹_Raman, is obtained from Ŵ^2−Ph. Hence, the Ŵ calculation of the entire relaxation time is τ⁻¹ = τ⁻¹_Orbach + τ⁻¹_Raman. For Raman, we have rotated the molecule in the Eigen frame of g-tensors of the ground state KD and then applied a small magnetic field to lift the degeneracy.

Results and discussion

Structure, bonding and the static relaxation mechanism in DyScS@C₈₂

To start, we optimized the geometry of DyScS@C₈₂ using DFT methods (see Computational Details). The optimized structure reveals a Dy–S bond distance of 2.514 Å (Fig. 1(a)), which falls within the range reported for the X-ray structure, despite its high degree of disorder, with a Dy–S–Sc angle of 119.8°. Additionally, we optimized the pentalene clusters DyScS(C₈H₆)₂ to assess the influence of the cage on the structure and magnetic properties.⁷⁰ In DyScS(C₈H₆)₂, the Dy–S bond distances is 2.584 Å, with a Dy–S–Sc angle of 120.4°, closely matching the value observed in the cage (Fig. 1(b)).

Fig. 1 (a) DFT optimised structure of DyScS@C₈₂, (b) DFT optimised structure of DyScS(C₈H₆)₂, (c) ab initio computed magnetic blockade diagram in DyScS@C₈₂. The thick black line indicates the KDs as a function of the computed magnetic moment. The green/blue arrows show the possible pathway through Orbach/Raman relaxation. The dotted red lines represent the presence of QTM/TA-QTM between the connecting pairs. The numbers provided at each arrow are the mean absolute values for the corresponding matrix element of the transition magnetic moment. Colour code: Dy^III -pink, S-yellow, C-grey and Sc-violet. Bond lengths and bond angle are given in Å and °, respectively.

To further investigate bonding characteristics, we performed Atoms in Molecules (AIM) analysis within the DFT framework (see Table S2† for details). In both DyScS@C₈₂ and DyScS(C₈H₆)₂, the Dy^III–cage and Dy^III–S interactions exhibit a |V(r)/G(r)| ratio of 1.15–1.28, suggesting a bonding nature that lies between covalent and ionic. Additionally, the positive Laplacian of the electron density (∇²ρ(r) > 0) suggests ionic interactions, whereas the negative total energy density (H(r) < 0) implies covalent character.⁷¹ Consequently, these interactions can be classified as polar covalent. To further evaluate metal–cage interactions, we calculated the delocalization index δ(Dy, C_cage), which falls within the range of 0.101–0.345. This is consistent with typical metal–cage interactions observed in EMFs (Fig. S1 and Table S3†).⁷²

Furthermore, we performed SA-CASSCF/RASSI/SINGLE_ANISO calculations using MOLCAS to evaluate the magnetic properties of DyScS@C₈₂ and DyScS(C₈H₆)₂. For DySc@C₈₂, the splitting of eight Kramers doublets (KDs) is found to be 1011.5 cm⁻¹. The g_xx, g_yy, and g_zz values of the ground KD are 0.002, 0.003, and 19.925 (Table S4†), respectively, indicating strong Ising anisotropy. The angle between the ground KD g_zz and the excited state g_zz remains small up to the second excited KD (2.5–5.9°) but increases to 20.9° in the third excited KD (Table S4†).

The third excited KD exhibits very strong transverse anisotropy (g_xx = 0.641, g_yy = 1.288, and g_zz = 8.531), leading to relaxation from this state with an estimated theoretical barrier for magnetization reversal (U_cal) of 800.4 cm⁻¹ (Fig. 1(c)). In contrast, the g tensors of DyScS(C₈H₆)₂ remain Ising-like up to the third excited state. Here, relaxation also occurs from the third excited state, but with a significantly lower energy barrier of 520.3 cm⁻¹ (Fig. S2†). This difference may be attributed to the shorter Dy–S bond distance in the cage structures, highlighting the crucial role of the cage in enhancing performance. Here it is noteworthy that the magnetic relaxation diagram computed using this approach is derived from transition matrix elements of the magnetic moment between opposite magnetization states (n+ → n−) provides and it provides only a qualitative picture of the magnetic relaxation.

Spin phonon relaxation in DyScS@C₈₂

For first-order transitions, the rate at which spin states |a〉 and |b〉 change is given by an eqn (5), first and second terms represent spin transitions caused by phonon absorption and emission, respectively. In perfectly axial Kramer systems without external time-reversal symmetry breaking interactions, direct J_z = 15/2 → J_z = −15/2 transitions are prohibited. Relaxation based on the given equation must occur through an excited state via phonon absorption. Whereas, second-order transitions involve two phonons simultaneously facilitating spin relaxation through an intermediate spin state |c〉 within the J = 15/2 manifold. The population transfer rate between spin states |a〉 and |b〉 in this case is described by another eqn (8). Here, G^2−Ph accounts for energy conservation and phonon thermal population in two-phonon interactions, similar to G^1−Ph for single-phonon processes. This process is mediated by the excited spin state |c〉. The spin-relaxation time is calculated using the lattice force constants, spin-phonon coupling coefficients (∂Ĥ_s/∂Q_α), in eqn (4) and (6), as mentioned earlier. Fig. 2 displays results from both first- and second-order perturbation theory, alongside the best fit to experimental data. First-order theory shows the expected exponential relationship between spin relaxation and temperature (T). Raman relaxation, however, the effective reversal barrier U_eff and the pre-exponential factor (τ₀) were determined by fitting the simulated Orbach data to the Arrhenius expressionIn contrast, the results for Raman relaxation adhere to a more complex mathematical relationship. Recent literature has suggested that Raman relaxation is expected to follow a specific temperature-dependent law.

Fig. 2 The relaxation times from the Orbach and Raman processes in comparison to the experimental relaxation times.

But if there is only a single pair of phonons are involved, where one phonon excites the system and another facilitates relaxation, the Raman relaxations are also simplified to

Based on Redfield equations, spin-phonon coupling coefficients are utilised to compute the spin-phonon relaxation time. One- and two-phonon processes have been simulated using second- and fourth-order density matrix time-dependent perturbation theory.^8,69,73 These simulations are conducted using the open-source programme MolForge, which can be downloaded from github.com/LunghiGroup/MolForge.

High temperatures see the dominance of the Orbach process, whereas the Raman process prevails at low temperatures. In Fig. 2, by plotting ln(τ) vs. 1/T, the slope of the curve at high temperatures gives the U_eff value of 264 cm⁻¹, which also coincides with the energy of the first excited Kramers Doublet (275.2 cm⁻¹). The U_cal value obtained from CASSCF calculations using the Molcas suite differs from this result because the previous method does not account for the effects of spin-phonon relaxation in its calculations. Although not identical, the similar molecule Dy₂S@C₈₂ studied by Chen, Krylov et al.^74,75 exhibits exchange-coupled behaviour at extremely low temperatures. However, at moderate and high temperatures, it demonstrates single-ion properties, making it suitable for comparison. Notably, it shows a similar high-temperature U_eff value of 363 cm⁻¹, which can be compared to the U_eff of 264 cm⁻¹ obtained using our method.

This aligns with the definition of the Orbach relaxation mechanism, which stipulates that phonons must be in resonance with the electronic states. However, the Orbach curve does not remain linear at low temperatures, indicating that other relaxation mechanisms, such as Raman relaxation, may become relevant. This is understandable because, at high temperatures, higher excited states are occupied, and phonons resonating with these states lead to relaxation. On the other hand, at low temperature a curve between ln(τ) vs. 1/T for Raman relaxation, the W_eff value comes out to be 49 cm⁻¹, which coincides with the first optical mode (45.7 cm⁻¹) at the Γ-point. Recent literature supports the idea that W_eff should align with the lowest-energy phonons strongly coupled to the magnetic moment.^{8,23,67,76,77} A similar phenomenon has been observed in metallofullerenes in the literature as well.^23,77 In previous papers, some of us observed that this phonon typically corresponds to one of the lowest energy modes.¹¹ This value is also comparable to the U_eff of 43 cm⁻¹ observed for Dy₂S@C₈₂ at moderate temperatures^74,75 and pointing to an interpreation of the latter as arising from Raman relaxation instead of Orbach as often assumed for these small barriers.¹¹ Importantly, our simulations show that for this compound, the crossover from Orbach to Raman relaxation happens at very high temperature and that Raman relaxation completely dominates the spin dynamics.

We note that the calculated relaxation times differ from the experimental results, except at high temperatures, as shown in Fig. 2. The use of the gas-phase approximation might partially explain the deviations between simulations and experiments at low temperature. Another important factor to consider is the omission of QTM contributions in our approach. Even though experiments have been conducted in the presence of a magnetic field, QTM might have not been completely quenched, making comparison to spin-phonon simulations not straightforward. Ideally a systematic experimental study as function of magnetic dilution could provide a more robust comparison. Despite these deviations at low temperature, experiments and simulations start converging at the highest temperatures at which experiments are available, pointing to intra-molecular contributions as the dominant in this regime.

To gain deeper insights into the relaxation mechanisms of the DyScS@C₈₂ complex, we measured transition probabilities for Orbach and the Raman relaxation, as shown in Fig. 3. Fig. 3(a) presents computed transition rates for the Orbach mechanism at 20 K, while Fig. 3(b) illustrates Raman relaxation probabilities at 6 K. It is noteworthy that these values represent the true ab initio computed spin-phonon transition rates used to calculate relaxation times, rather than the typical expected dipole moment values associated with such mechanisms. These diagrams suggest that Orbach mechanism is promoted through the absorption of a phonon resonant with the first excited KD, while Raman relaxation is instead described by an intra-KD transition promoted by the simultaneous absorption and emission of pairs of degenerate phonons.

Fig. 3 (a) Orbach transition rates computed at 20 K and (b) Raman transition rates computed at 6 K. The Orbach transition rates show the matrix elements of, Ŵ^1−ph whereas the Raman relaxations show the matrix elements of Ŵ^2−ph. Rates are expressed in ps⁻¹ and only rates larger than 10⁻¹⁶ are reported explicitly. The x-axis displays the computed average magnetic moment for the first four KDs and their energy separation from the ground state.

Examination of vibrational density of states, spin-phonon coupling constants and the molecular vibrations

After establishing the temperature dependence of the relaxation times, we focused on the spin-phonon coupling constants and the vibrational density of states (DOS) to decipher the origin of the Raman relaxation observed in the aforementioned section. Fig. 4(a) illustrates the spin-phonon coupling coefficients as a function of frequencies along with the vibrational DOS. Fig. 4(b) presents the detailed spin-phonon coupling coefficients alongside the energies of the eight Kramers doublets (KDs). We found no direct correlation between the vibrational DOS and spin-phonon coupling coefficients at low frequencies, though a partial correlation is evident at higher frequencies. Two regions exhibit very strong spin-phonon coupling: one near 50 cm⁻¹ (vib1, shown in Fig. 5(a and b)) and another at 87 cm⁻¹ (vib2, shown in Fig. 5(c and d)). Both modes occur at significantly lower frequencies than the first excited KD, indicating no possibility of Orbach relaxation for this molecule. Closer inspection of these modes revealed that they involve the movement of the Dy–S–Sc moiety inside the cage. Other spin-phonon coupling coefficients at higher frequencies are much smaller.

Fig. 4 (a) Spin-phonon coupling distribution and the vibrational density of states as functions of energy. A Gaussian smearing with σ = 1 cm⁻¹ and σ = 10 cm⁻¹ has been applied to the two functions, respectively. (b) Close-up of the spin-phonon coupling distribution overlapped with the energies of the KDs.

Fig. 5 (a) First vibration of vib1, (b) second vibration of vib1 (c) first vibration of vib2 (d) second vibration of vib2. The red arrows indicate the movement, with the size of the arrow representing the extent of motion for that particular atom.

As expected, since anisotropy should be primarily affected by the first coordination sphere, the highest spin-phonon coupling is observed with the movements of S, Sc, and Dy atoms inside the fullerene. We will now discuss vib1 and vib2 in detail. Both modes are labelled in Fig. 4(a) and are pictorially represented in Fig. 5. These modes were selected based on their highest spin-phonon coupling coefficients and exhibit similar types of motions. vib1 consists of the first two phonons at 45.7 cm⁻¹ and 54.0 cm⁻¹, involving two different vibrations. This aligns with our W_eff value (49 cm⁻¹), as well as the low temperature energy barrier of the Dy₂S@C₈₂. First is a wagging-like vibration, where the Dy–S−Sc angle (and Sc–S, Dy–S bond distances) changes from 119.8° (2.331 Å, 2.514 Å) to 115.8° (2.395 Å, 2.550 Å). The second vibration is more like a scissoring motion, with Dy^III moving more, changing the angle from 115.4° (2.333 Å, 2.532 Å) to 124.4° (2.328 Å, 2.518 Å). The vib2 consists of two phonons at 79.8 cm⁻¹ and 90.4 cm⁻¹. The first phonon corresponds to a similar wagging mode as in vib1, where the bond angle and bond distances change from 117.7° (2.357 Å, 2.536 Å) to 119.8° (2.330 Å, 2.514 Å). The second phonon in vib2 corresponds to a scissoring-like bending motion, with Sc^III moving more, changing the angle from 117.1° (2.347 Å, 2.509 Å) to 122.7° (2.335 Å, 2.524 Å).

Additionally, it was observed that cage atoms begin to visibly vibrate only at the eighth phonon mode, which occurs at 210.9 cm⁻¹. This corresponds to the previously discussed U_eff value (264 cm⁻¹) and the one reported for Dy₂S@C₈₂ in the context of Orbach relaxation. This suggests that the cage indeed provides a largely rigid ligand environment that does not contribute to Raman relaxation.

Discussion and conclusions

The field of molecular magnetism has made significant progress over the past three decades, with the discovery of compounds exhibiting remarkably U_eff values and magnetic hysteresis at increasingly high temperatures. However, the field now faces a critical juncture, as traditional strategies for increasing the effective energy barrier in single-ion complexes may be approaching their practical limits. Future advancements in designing molecular compounds with extended spin lifetimes will likely stem from a more comprehensive understanding of the entire spin relaxation process. This necessitates a shift in focus beyond the singular pursuit of higher U_eff values, encompassing a broader analysis of relaxation phenomena, including the often overlooked preexponential factor τ₀. On the other hand, first-principles simulations have recently demonstrated their capacity to accurately model spin relaxation processes, enabling the estimation of spin relaxation times without necessarily relying on experimental data and providing an unbiased means to benchmark our understanding of spin-phonon interactions and relaxation mechanisms. This computational approach enables researchers to directly interpret relaxation experiments without relying on phenomenological models. By providing a more fundamental and system-specific understanding, these simulations offer a powerful alternative to traditional interpretative frameworks, potentially leading to more accurate and insightful analyses of spin relaxation phenomena in molecular magnetic materials.

In this paper, we have delved into the spin-phonon relaxation mechanisms within fullerene molecules, focusing on the high-performance SMM DyScS@C₈₂, which shows a T_B of 7.3 K. The study of spin-phonon relaxation in DyScS@C₈₂ SMMs offers significant insights into how the interactions between molecular vibrations and electronic spin states influence the magnetic behaviour of these systems. The research highlights the importance of a molecular-level understanding of blocking temperatures in SMMs, which has long been a critical challenge in advancing practical applications of these molecules. The key-findings of this work are outlined below.

(i) Role of cage in SMM performance

The fullerene cage plays a crucial role in stabilizing the molecular fragment and modifying its geometry to enhance SMM performance. The spatial constraints imposed by the cage shorten Dy-C_cage and Dy–S bond distances, which, while not enhancing axiality through the ligand field, lead to a higher effective energy barrier (U_cal). This effect is evident in DyScS@C₈₂ (800.4 cm⁻¹) compared to DyScS(C₈H₆)₂ (520.3 cm⁻¹), demonstrating that the cage-induced bond contraction contributes significantly to improved SMM behavior.

(ii) Mode of relaxation

Analysis of relaxation time versus temperature graphs reveals that the predominant relaxation mechanism is the Raman process. The Orbach mechanism involves phonon absorption resonant with the first excited KD, while Raman relaxation occurs via intra-KD transitions through simultaneous absorption and emission of degenerate phonons. The highest spin-phonon coupling coefficients occur at low frequencies, significantly lower than the energy of the first Kramers doublet (KD), suggesting that the relaxation mechanism cannot be attributed to the Orbach process.

(iii) Spin-phonon coupling dynamics in DyScS@C₈₂

At lower frequencies, we observed a pronounced phenomenon of spin-phonon coupling, where interactions between spin states and vibrational modes are notably strong. The lowest energy optical phonons arise from specific motions within the {Dy–S−Sc} moiety, exhibiting a robust coupling with the spin states of the system. It is noteworthy that vibrations originating from the fullerene cage occur at higher frequencies, beginning prominently at 210.9 cm⁻¹. These vibrations primarily represent the structural dynamics of the cage itself, with minimal interaction with the enclosed Dy–S−Sc moiety, thereby not leading to significant spin-phonon coupling. Thus, the cage clearly offers a significant advantage and does not contribute to spin-phonon relaxation.

(iv) How cluster atom vibrations drive relaxation

We have identified four distinct low-lying phonons crucial to the relaxation dynamics of the system. These phonons involve specific bending motions that cause variations in the Dy–S–Sc bond angles, thereby initiating relaxation processes. The changes in Dy–S–Sc bond angles are known to influence the axial ligand field: an increase in bond angle reduces the axial ligand field strength, consequently promoting transverse anisotropy and facilitating relaxation pathways within the system. Unfortunately, these phonons predominantly operate at low frequencies, which poses a limitation on the system's ability to achieve high blocking temperatures. This characteristic underscores the importance of understanding and managing these low-frequency vibrational modes in the design and engineering of materials aimed at enhancing thermal stability and magnetic properties.

Our study reveals that while the fullerene cage improves stability and crystal field splitting of the {DyScS} fragment, it does not completely eliminate low-frequency vibrational modes within the cluster, which contribute to spin-phonon relaxation. These key low-lying phonons affect Dy–S–Sc bond angles, leading to transverse anisotropy and relaxation pathways. For enhanced performance in future EMFs, designs should focus on stronger Dy-ligand bonds, smaller cage sizes to reduce internal tumbling, and maintaining the cluster close to its equilibrium geometry. This is exemplified in the literature evidences for instance, DyLu₂N@I_h(7)-C₈₀, with a blocking temperature of 9 K,³⁰ demonstrates stronger Dy–ligand interaction, while optimal combinations like in Dy@Cs(6)-C₈₁N,²⁶ show fewer low-energy phonons coupled to spin states, leading to even higher blocking temperatures. This novel approach not only clarifies the spin-phonon mechanism for EMF-based SMMs but also provides a strategic pathway to enhance the performance of future EMF-based SMMs. Overall, this study advances the understanding of spin-phonon relaxation in SMMs and provides valuable guidelines for future molecular design aimed at increasing T_B values. By focusing on the interplay between molecular structure, vibrational modes, and spin dynamics, it opens up new avenues for the development of high-performance SMMs, particularly in the context of fullerene-based systems where rigid molecular environments can stabilize low-coordinate lanthanide centers.

Data availability

The data supporting this article have been included as part of the ESI.†

Author contributions

GR and AL designed and supervised the investigations; TS carried out DFT, periodic DFT, CASSCF, and ab initio spin-phonon relaxation calculations; RKT and SD contributed to the CASSCF calculations; LAM assisted in troubleshooting the ab initio spin-phonon relaxation calculations.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

TS is thankful to CSIR India for funding. GR would like to thank SERB (SB/SJF/2019-20/12; CRG/2022/001697) and SCPP (RD/0523-IOE0010-007) for funding. AL is thankful to the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 948493).

References

1. F.-S. Guo, B. M. Day, Y.-C. Chen, M.-L. Tong, A. Mansikkamäki and R. A. Layfield, Science, 2018, 362, 1400–1403.
2. C. A. Goodwin, F. Ortu, D. Reta, N. F. Chilton and D. P. Mills, Nature, 2017, 548, 439–442.
3. S. K. Gupta, T. Rajeshkumar, G. Rajaraman and R. Murugavel, Chem. Sci., 2016, 7, 5181–5191.
4. A. B. Canaj, S. Dey, E. R. Marti, C. Wilson, G. Rajaraman and M. Murrie, Angew. Chem., Int. Ed., 2019, 58, 14146–14151.
5. T. Sharma, M. K. Singh, R. Gupta, M. Khatua and G. Rajaraman, Chem. Sci., 2021, 12, 11506–11514.
6. C. A. Gould, K. R. McClain, D. Reta, J. G. Kragskow, D. A. Marchiori, E. Lachman, E.-S. Choi, J. G. Analytis, R. D. Britt and N. F. Chilton, Science, 2022, 375, 198–202.
7. A. Lunghi, F. Totti, R. Sessoli and S. Sanvito, Nat. Commun., 2017, 8, 14620.
8. M. Briganti, F. Santanni, L. Tesi, F. Totti, R. Sessoli and A. Lunghi, J. Am. Chem. Soc., 2021, 143, 13633–13645.
9. J. G. Kragskow, A. Mattioni, J. K. Staab, D. Reta, J. M. Skelton and N. F. Chilton, Chem. Soc. Rev., 2023, 4567–4585.
10. M. Dolg, U. Wedig, H. Stoll and H. Preuss, J. Chem. Phys., 1987, 86, 866–872.
11. S. Mondal and A. Lunghi, J. Am. Chem. Soc., 2022, 144, 22965–22975.
12. D. N. Woodruff, R. E. Winpenny and R. A. Layfield, Chem. Rev., 2013, 113, 5110–5148.
13. L. Ungur and L. F. Chibotaru, Inorg. Chem., 2016, 55, 10043–10056.
14. L. Ungur and L. F. Chibotaru, Phys. Chem. Chem. Phys., 2011, 13, 20086–20090.
15. J. D. Rinehart and J. R. Long, Chem. Sci., 2011, 2, 2078–2085.
16. A. Ullah, J. Cerda, J. J. Baldovi, S. A. Varganov, J. Arago and A. Gaita-Arino, J. Phys. Chem. Lett., 2019, 10, 7678–7683.
17. A. Lunghi and S. Sanvito, Sci. Adv., 2019, 5, eaax7163.
18. T. Sharma, R. Rana, A. Swain and G. Rajaraman, Inorg. Chem. Front., 2025 10.1039/d5qi00907c.
19. J. P. Durrant, J. Tang, A. Mansikkamäki and R. A. Layfield, Chem. Commun., 2020, 56, 4708–4711.
20. R. Nabi, R. K. Tiwari and G. Rajaraman, Chem. Commun., 2021, 57, 11350–11353.
21. G. Velkos, D. Krylov, K. Kirkpatrick, X. Liu, L. Spree, A. Wolter, B. Büchner, H. Dorn and A. Popov, Chem. Commun., 2018, 54, 2902–2905.
22. M. Nie, J. Xiong, C. Zhao, H. Meng, K. Zhang, Y. Han, J. Li, B. Wang, L. Feng and C. Wang, Nano Res., 2019, 12, 1727–1731.
23. L. Spree, C. Schlesier, A. Kostanyan, R. Westerström, T. Greber, B. Büchner, S. M. Avdoshenko and A. A. Popov, Chem. –Eur. J., 2020, 26, 2436–2449.
24. Y. Wang, J. Xiong, J. Su, Z. Hu, F. Ma, R. Sun, X. Tan, H.-L. Sun, B.-W. Wang and Z. Shi, Nanoscale, 2020, 12, 11130–11135.
25. G. Velkos, W. Yang, Y.-R. Yao, S. M. Sudarkova, F. Liu, S. M. Avdoshenko, N. Chen and A. A. Popov, Chem. Commun., 2022, 58, 7164–7167.
26. Z. Hu, Y. Wang, A. Ullah, G. M. Gutiérrez-Finol, A. Bedoya-Pinto, P. Gargiani, D. Shi, S. Yang, Z. Shi and A. Gaita-Ariño, Chem, 2023, 9, 3613–3622.
27. W. Yang, G. Velkos, M. Rosenkranz, S. Schiemenz, F. Liu and A. A. Popov, Adv. Sci., 2024, 11, 2305190.
28. C. Schlesier, L. Spree, A. Kostanyan, R. Westerström, A. Brandenburg, A. U. Wolter, S. Yang, T. Greber and A. A. Popov, Chem. Commun., 2018, 54, 9730–9733.
29. R. Westerstrom, J. Dreiser, C. Piamonteze, M. Muntwiler, S. Weyeneth, H. Brune, S. Rusponi, F. Nolting, A. Popov and S. Yang, J. Am. Chem. Soc., 2012, 134, 9840–9843.
30. Y. Hao, G. Velkos, S. Schiemenz, M. Rosenkranz, Y. Wang, B. Büchner, S. M. Avdoshenko, A. A. Popov and F. Liu, Inorg. Chem. Front., 2023, 10, 468–484.
31. D. Krylov, F. Liu, S. Avdoshenko, L. Spree, B. Weise, A. Waske, A. Wolter, B. Büchner and A. Popov, Chem. Commun., 2017, 53, 7901–7904.
32. A. Kostanyan, C. Schlesier, R. Westerström, J. Dreiser, F. Fritz, B. Büchner, A. A. Popov, C. Piamonteze and T. Greber, Phys. Rev. B, 2021, 103, 014404.
33. G. Velkos, W. Yang, Y.-R. Yao, S. M. Sudarkova, X. Liu, B. Büchner, S. M. Avdoshenko, N. Chen and A. A. Popov, Chem. Sci., 2020, 11, 4766–4772.
34. W. Yang, G. Velkos, F. Liu, S. M. Sudarkova, Y. Wang, J. Zhuang, H. Zhang, X. Li, X. Zhang and B. Büchner, Adv. Sci., 2019, 6, 1901352.
35. W. Yang, G. Velkos, S. Sudarkova, B. Büchner, S. M. Avdoshenko, F. Liu, A. A. Popov and N. Chen, Inorg. Chem. Front., 2022, 9, 5805–5819.
36. A. Brandenburg, D. S. Krylov, A. Beger, A. U. Wolter, B. Büchner and A. A. Popov, Chem. Commun., 2018, 54, 10683–10686.
37. W. Cai, J. D. Bocarsly, A. Gomez, R. J. L. Lee, A. Metta-Magaña, R. Seshadri and L. Echegoyen, Chem. Sci., 2020, 11, 13129–13136.
38. B. G. Lippert, J. H. Parrinello and P. Michele, Mol. Phys., 1997, 92, 477–488.
39. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
40. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt and F. Schiffmann, J. Chem. Phys., 2020, 152.
41. G. Lippert, J. Hutter and M. Parrinello, Theor. Chem. Acc., 1999, 103, 124–140.
42. J. Paier, R. Hirschl, M. Marsman and G. Kresse, J. Chem. Phys., 2005, 122, 234102.
43. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127.
44. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
45. D. A. Pantazis and F. Neese, J. Chem. Theor. Comput., 2009, 5, 2229–2238.
46. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
47. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
48. P. A. Malmqvist, B. O. Roos and B. Schimmelpfennig, Chem. Phys. Lett., 2002, 357, 230–240.
49. F. Aquilante, J. Autschbach, R. K. Carlson, L. F. Chibotaru, M. G. Delcey, L. De Vico, I. F. Galvan, N. Ferre, L. M. Frutos and L. Gagliardi, J. Comput. Chem., 2016, 37, 506–541.
50. L. F. Chibotaru and L. Ungur, J. Chem. Phys., 2012, 137.
51. B. r. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov, P.-O. Widmark and A. C. Borin, T J. Chem. Phys. A, 2008, 112, 11431–11435.
52. S. Dey, T. Sharma and G. Rajaraman, Chem. Sci., 2024, 15, 6465–6477.
53. S. Dey and G. Rajaraman, Chem. Sci., 2021, 12, 14207–14216.
54. M. Frisch, G. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. Petersson, Gaussian 09, Revision d. 01, Gaussian Inc., Wallingford CT, 2009, vol. 201.
55. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
56. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785.
57. A. D. Becke, Phys. Rev. A, 1988, 38, 3098.
58. D. Andrae, U. Haeussermann, M. Dolg, H. Stoll and H. Preuss, Theor. Chim. Acta, 1990, 77, 123–141.
59. T. R. Cundari and W. J. Stevens, J. Chem. Phys., 1993, 98, 5555–5565.
60. M. Coletta, S. Sanz, D. J. Cutler, S. J. Teat, K. J. Gagnon, M. K. Singh, E. K. Brechin and S. J. Dalgarno, Dalton Trans., 2020, 49, 14790–14797.
61. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
62. S. Haldar, L. A. Mariano, A. Lunghi and L. Gagliardi, J. Chem. Theor. Comput., 2025, 21, 2829–2838.
63. A. Abragam and B. Bleaney, Electron paramagnetic resonance of transition ions, OUP Oxford, 2012.
64. L. Ungur and L. F. Chibotaru, Chem. - Eur. J., 2017, 23, 3708–3718.
65. D. Reta, J. G. Kragskow and N. F. Chilton, J. Am. Chem. Soc., 2021, 143, 5943–5950.
66. K.-X. Yu, J. G. Kragskow, Y.-S. Ding, Y.-Q. Zhai, D. Reta, N. F. Chilton and Y.-Z. Zheng, Chem, 2020, 6, 1777–1793.
67. A. Lunghi and S. Sanvito, J. Chem. Phys., 2020, 153.
68. R. Nabi, J. K. Staab, A. Mattioni, J. G. Kragskow, D. Reta, J. M. Skelton and N. F. Chilton, J. Am. Chem. Soc., 2023, 145, 24558–24567.
69. A. Lunghi, Sci. Adv., 2022, 8, eabn7880.
70. D. Qingming, J. Am. Chem. Soc., 2014, 136, 4257–4264.
71. P. Macchi and A. Sironi, Coord. Chem. Rev., 2003, 238, 383–412.
72. A. A. Popov and L. Dunsch, Chem. –Eur. J., 2009, 15, 9707–9729.
73. A. Lunghi and G. Rajaraman, Challenges Adv. Comput. Chem. Phys., 2023, 219–289.
74. C.-H. Chen, D. S. Krylov, S. M. Avdoshenko, F. Liu, L. Spree, R. Yadav, A. Alvertis, L. Hozoi, K. Nenkov and A. Kostanyan, Chem. Sci., 2017, 8, 6451–6465.
75. D. Krylov, G. Velkos, C.-H. Chen, B. Büchner, A. Kostanyan, T. Greber, S. M. Avdoshenko and A. A. Popov, Inorg. Chem. Front., 2020, 7, 3521–3532.
76. L. Gu and R. Wu, Phys. Rev. Lett., 2020, 125, 117203.
77. Z. Hu, A. Ullah, H. Prima-Garcia, S. H. Chin, Y. Wang, J. Aragó, Z. Shi, A. Gaita-Ariño and E. Coronado, Chem. –Eur. J., 2021, 27, 13242–13248.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sc07786e

---