Introduction

Lanthanides are well-known for their large atomic moments and magnetic anisotropies, and embedding discrete lanthanide ions in molecular environments leads to nanomagnets exhibiting magnetic bistability and slow relaxation of magnetization on a single molecule level1,2. Single molecule magnets (SMMs) can be used as core elements of nanospintronic devices3, such as spin valves4, spin transistors5,6, or building blocks of quantum computers7,8. Optimization of symmetry and ligand environment led to a dramatic improvement of lanthanide SMMs during the last decade9,10,11 with the latest discovery of magnetic hysteresis at 60–80 K in Dy-metallocenium salts12,13,14,15. Effective barriers of magnetization reversal higher than 1000 K were reported for several lanthanide SMMs12,13,14,15,16,17,18,19,20. Single lanthanide atoms were also found to keep magnetic bistability up to 45 K on MgO|Ag(100) substrate21,22,23.

Combination of several magnetic centers within one molecule may lead to high-spin ground states and can largely suppress quantum tunneling, which is the main low-temperature relaxation mechanism for single-ion magnets in zero magnetic field. Design of multinuclear SMMs has been a viable strategy since the discovery of SMM behavior in the Mn12 complex24. The temperature range at which a multinuclear magnet can be considered as a giant spin rather than a combination of weakly interacting individual spins is limited by the strength of exchange interactions. Whereas exchange interactions between transition metals can be tuned in a wide range, the localized nature of 4f electrons results in weak exchange interactions in lanthanide compounds rarely exceeding 1 cm−1. As a result, when the relaxation of magnetization in multinuclear lanthanide SMMs is driven via exchange excitations, the barriers to magnetization reversal are usually well below 100 K25.

Exchange coupling in lanthanide molecule magnets can be increased by introducing radical bridges26. The radical bridge usually features a rather diffuse singly occupied molecular orbital, which exhibits stronger interactions with the 4f electrons. The lanthanide-radical exchange coupling constants can reach values of −27 cm−1. The strongest coupling so far has been found in dilanthanide complexes with N23− radical bridges27,28,29, and the corresponding Tb complex has a blocking temperature near 30 K, which is the highest known value among multinuclear SMMs30. In this work we explore the limiting case of this concept in which the role of a radical bridge is played by a single unpaired electron, residing on the lanthanide-lanthanide bonding orbital and coupling the lanthanide spins inside a fullerene.

The empty space inside carbon cages provides unmatched possibilities for stabilizing small metallic clusters in unconventional valence and spin states31, such as metal dimers in dimetallofullerenes. Although Coulomb repulsion between two metal ions prevails over the covalent bonding32,33,34, the metal dimers cannot dissociate inside fullerenes. This unique situation allows for direct Ln–Ln bonding in dimetallofullerenes35,36,37,38, which could not be realized in any other molecular lanthanide compound so far. Of particular interest are Ln2@C80 molecules. The C80-Ih fullerene cage is unstable in the neutral state due to the presence of a four-fold degenerate orbital occupied by only two electrons. At the same time, the hexaanion of C80-Ih has a very stable closed-shell electronic structure, and hence this fullerene is the most preferable host for endohedral species acting as donors of six electrons, such as the early lanthanide dimers La2 or Ce239. However, in the Ln2@C80 molecules with heavier lanthanides (Gd to Lu), the Ln2 dimers transfer only 5 electrons to the fullerene orbitals, leaving one electron on the Ln–Ln bonding orbital40. The electronic structure of the fullerene cage in such Ln25+@C805− dimetallofullerenes can be stabilized by addition of one electron41, substitution of one carbon by a nitrogen atom (giving azafullerenes Ln2@C79N37,42) or by functionalization with a radical group (giving derivatives Ln2@C80R, R = CF340,43, CH2Ph38). A formal oxidation state of the lanthanides in such dimetallofullerenes is Ln+2.5. EPR studies of Gd2@C79N revealed a ferromagnetic ground state with a giant spin of S = 15/2, and magnetization studies showed that the exchange coupling between the Gd spin and the unpaired electron residing on the Gd−Gd bonding orbital is as large as 170 cm−1 (refs. 44,45). Very recently we have reported that Dy2@C80(CH2Ph) shows exceptional magnetic properties with magnetic hysteresis up to 22 K and Dy-electron exchange coupling of 32 cm−1 38.

In order to investigate and understand the principles underlying the SMM behavior in encapsulated lanthanide dimers, in this work we synthesize and study an array of Ln2@C80(CH2Ph) molecules (Ln2 = Y2, Gd2, Tb2, Dy2, Ho2, Er2, TbY, TbGd), all featuring single-electron Ln–Ln bonding molecular orbitals (MO). This bonding situation leads to giant exchange interactions in all magnetic molecules and is very beneficial for the molecular magnetism, especially in Tb2@C80(CH2Ph) showing a giant coercivity and the highest known blocking temperature among dinuclear lanthanide complexes. Furthermore, we demonstrate that the single-electron Ln–Ln bonding MO is redox-active, which allows control of the magnetism via electron transfer.

Results

Synthesis and internal dynamics

All compounds were synthesized using the extraction/functionalization procedure recently developed in our group (Supplementary Fig. 1)38. Carbon soot obtained by arc-discharge of metal-oxide filled graphite rods is extracted with hot DMF, giving a mixture of fullerene anions in DMF solution. The anions are reacted with benzyl-bromide to yield neutral air-stable benzyl monoadducts, which are further separated by multistep HPLC protocol to yield pure Ln2@C80(CH2Ph) derivatives (denoted as {Ln2} hereafter, Ln = Y, Gd, Tb, Dy, Ho, Er; Supplementary Figs 220). If a mixture of two metal oxides (Ln and Ln′) is used in the arc-discharge synthesis, the same procedure gives a mixture of {Ln2}, {Ln′2}, and {LnLn′}. For the Tb-Y mixed-metal system, the mixture is further separated by recycling HPLC to give pure {Y2}, {Tb2} and {TbY}. In the Tb-Gd system, the separation into individual components with HPLC could not be achieved, and the studied sample comprised ca 20% {Gd2}, 30% {Tb2} and 50% {TbGd}. Despite the unconventional oxidation state of the lanthanides (+2.5) and the presence of an unpaired valence electron (Fig. 1), the {Ln2} compounds are air-stable at room temperature (Supplementary Fig. 21) and do not require special handling conditions as many Ln compounds in unconventional oxidation states do.

Fig. 1
figure 1

Molecular structure of Ln2@C80(CH2Ph). Single-occupied Ln–Ln bonding molecular orbital (left; carbons are gray, hydrogens are white, lanthanides are green), and schematic depiction of the molecule (right; the arrow indicates an unpaired electron residing on the Ln–Ln bonding orbital)

Molecular structure and internal dynamics of the metal atoms inside the fullerene cage are studied by variable-temperature single-crystal X-ray diffraction for {Dy2} as a representative example providing single crystals of sufficient quality (Supplementary Figs 22, 23, Supplementary Table 1). Virtually identical Vis-NIR, IR, and Raman spectra prove that all isolated {Ln2} compounds are isostructural (Supplementary Figs 2426). Figure 2a shows the molecular structure of {Dy2} at 100 and 290 K along with the temperature dependence of atomic displacement parameters. Whereas the fullerene cage with the attached benzyl group remains ordered at all temperatures, the metal atoms show a pronounced increase of their mobility with temperature. At 100 K each Dy atom occupies two well-defined positions with 70 and 30% occupancy. The increase of temperature to 290 K induces vigorous motions of the metal atoms as can be seen from the increase of the number of metal sites and increase of their displacement parameters compared to those of the carbon atoms (Fig. 2a, Supplementary Table 2, Supplementary Note 1). The disordered metal positions tend to be distributed within a single plain, which suggests that the Ln2 dimer in {Ln2} exhibits in-plain rotation near room temperature.

Fig. 2
figure 2

Structure and dynamics of Ln2@C80(CH2Ph). a Molecular structure of {Dy2} at 100 and 290 K and atomic displacement parameters as a function of temperature between 100 and 290 K (to guide the eye, vertical lines separate displacement parameters of the C80 cage, CH2Ph group, and Dy atoms); Dy atoms are shown as spheres with radii proportional to the site occupancies; b Raman spectra of {Ln2} compounds in the low-frequency range measured at 77 K; two metal-cage stretching modes are indicated by red dotted lines, the most prominent fullerene cage squashing mode is indicated by a black dotted line, atomic displacements for the metal-cage stretching modes near 150 and 165 cm−1 are also shown on the right; c 1H NMR spectra of {Ln2} measured in CS2 solution at room temperature; {Tb2} signals in the spectrum of the {TbGd} sample are denoted by asterisks

Translational motions of the metal atoms inside the cage manifest in low-frequency vibrations in the Raman spectra (Fig. 2b). A strong Raman line near 150 cm−1 is due to the in-phase cage-metal stretching mode in {Ln2}. Its frequency exhibits pronounced metal dependence and shifts from 144.5 cm−1 in {Er2} to 151.2 cm−1 in {Gd2} and further to 181.1 cm−1 in {Y2} in accordance with the decrease of the metal atomic mass. Similar metal dependence is found for the anti-symmetric metal-cage mode observed as a low-intensity peak shifting from 163.7 cm−1 in {Er2} to 169.6 cm−1 in {Gd2}. At higher frequencies only carbon atoms contribute to the vibrational displacements, and the IR and Raman spectra of all studied {Ln2} compounds are virtually identical (Supplementary Figs 25, 26).

Paramagnetic 1H NMR spectroscopy

The Ln2 dimers inside the fullerene cages act as magnets, creating dipolar magnetic fields. The strength and spatial distribution of these fields can be evaluated by 1H NMR spectroscopy using benzyl protons as a probe. In solution 1H NMR spectra, all {Ln2} compounds except for {Y2} and {Gd2} exhibit well-defined temperature-dependent 1H resonances (Fig. 2c, Supplementary Figs 27, 33) strongly shifted from the standard chemical shifts of the benzyl group in diamagnetic compounds (3–7 ppm). These paramagnetic shifts are caused by the screening of the external magnetic field by the dipolar field of the endohedral lanthanide dimer. Since the molecules in solution are rotating fast, the isotropic contributions average out and the paramagnetic shift (δpara) serves as a measure of a magnetic anisotropy, taking the following form in a point dipole approximation:46

$$\delta _i^{{\mathrm{para}}} = \frac{{\left( {3{\mathrm{cos}}^2\theta _i - 1} \right)}}{{12{\mathrm{\pi }}R_i^3}}\left( {\chi _\parallel ^{{\mathrm{Ln}}_2} - \chi _ \bot ^{{\mathrm{Ln}}_2}} \right)$$
(1)

where the first term is defined via polar coordinates Ri and θi of the i-th proton in the coordinate system centered on the Ln2 dimer with polar axis along the Ln–Ln bond and is expected to be very similar for all {Ln2} compounds. The second term in Eq. 1 is the difference of the longitudinal \(\chi _\parallel\) and transverse \(\chi _ \bot\) magnetic susceptibilities of the [Ln3+e–Ln3+] system. {Dy2} and {Tb2} exhibit almost identical 1H chemical shifts, indicating a similarity of their magnetic properties, whereas the values of {Ho2} are at least twice smaller. Substitution of one Tb in {Tb2} by Y results in a 4.3-fold drop of the paramagnetic shift in {TbY}. Remarkably, 1H resonances in {Er2} are shifted into the positive direction, revealing that the sign of the magnetic anisotropy in the endohedral Er2 dimer is opposite to that in Tb2, Dy2, and Ho2 dimers.

The isotropic spin of Gd does not induce dipolar shifts in the NMR spectra (although it does affect the relaxation rates of proton spins and thus makes the lines very broad—we, therefore, could not detect the 1H NMR spectrum of {Gd2}). If the Gd spin in {TbGd} were behaving isotropically, the 1H chemical shifts of {TbGd} would be close to those of {TbY}. However, the measured 1H shifts in {TbGd} are two times larger than in {TbY} (Fig. 2c). This shows that the Gd spin is locked to the anisotropic Tb spin by exchange interaction through the unpaired electron spin. Observation of this phenomenon at room temperature indicates that the exchange coupling is very strong.

Coupling of lanthanide spins in {Ln2}

The combination of large anisotropic Ln spins with a strong exchange coupling via a delocalized unpaired electron in {Ln2} molecules is promising for SMM. The study of an array of {Ln2} molecules with different lanthanides enables disentanglement of the exchange and anisotropy factors in determining magnetic properties of the {Ln2} system. This task, however, requires a common theoretical framework, which is outlined in this section.

The effective spin Hamiltonian of the {Ln2} molecule includes single-ion ligand-field (LF) effects together with the Kondo description of magnetic interactions between the lowest J-multiplet of lanthanides with the partially delocalized unpaired electron occupying the spd-hybrid Ln–Ln bonding orbital:

$$\hat H_{{\mathrm{spin}}} = \hat H_{{\mathrm{LF}}_1} + \hat H_{{\mathrm{LF}}_2} + \hat H_{{\mathrm{sf}}},$$
(2)

where \(\hat H_{{\mathrm{LF}}_i}\) is the single-ion LF Hamiltonian of the i-th lanthanide site. Later, the LF Hamiltonian will be limited to the crystal-field shape commonly used for 4f systems. Further, \(\hat H_{{\mathrm{sf}}}\) is the spin-fermion Hamiltonian describing direct exchange interactions between the lanthanide ion moments, the kinetic energy for the electron, and an on-site exchange interaction:

$$\hat H_{{\mathrm{sf}}} = - 2j_{12}\hat J_{{\mathrm{Ln}}_1}\hat J_{{\mathrm{Ln}}_2} + t\mathop {\sum }\limits_\sigma \left[ {c_{1\sigma }^\dagger c_{2\sigma }^{} + c_{2\sigma }^\dagger c_{1\sigma }^{}} \right] - 2\hat s(K_1\hat J_{{\mathrm{Ln}}_1} + K_2\hat J_{{\mathrm{Ln}}_2}).$$
(3)

Here j12 is the direct exchange coupling between the localized lanthanide moments \(\hat J_{{\mathrm{Ln}}_i}\), t is the electron hopping amplitude between sites 1 and 2, c (\(c_{i\sigma }^\dagger\)) is the creation (annihilation) operator for the electron at site i with spin σ, and Ki denotes the on-site Kondo exchange coupling constant between the localized 4f moment \(\hat J_{{\mathrm{Ln}}_i}\) and the delocalized spin \(\hat s\). In the modelling of magnetic properties of {Ln2} molecules using spin Hamiltonians derived from Eqs. 2 and 3, lanthanide moments \(\hat J_{{\mathrm{Ln}}_i}\) are treated in a full \(|J,m_J\rangle\) basis set for each lanthanide ion. In general, spin-spin interactions in Eq. 3 are anisotropic and may lead to a rich phase diagram in the full parameter space due to an interplay of very different energy scales47,48. A complete derivation of all the parameters is beyond the scope of the present report49,50. Nevertheless, this spin Hamiltonian can be simplified by considering the symmetry and properties of individual {Ln2} systems.

Magnetic properties of {Gd2}

Magnetization curves of {Gd2} show no hysteresis down to 1.8 K, and we did not succeed in determining relaxation of magnetization by using AC magnetometry. Thus, {Gd2} is not a SMM, which is not very surprising since Gd3+ ions are magnetically isotropic. Hence, {Gd2} provides a convenient example to study the role of the exchange interactions without the contribution of any single-ion anisotropy.

A recent theoretical study of the [Gd2@C80] anion showed that the hopping amplitude t exceeds 10,000 cm−1 (ref. 51), which suggests a weak Kondo coupling limit (\(K/t \ll 1\)). At this limit, the kinetic term in Eq. 3 can be omitted, and the form of the spin Hamiltonian can be simplified to a simple 3-center model with only exchange interactions between the lanthanides and the electron spin30,44,45,49,52,53,54. As the ligand-field terms for isotropic Gd spins can be neglected as well, the effective spin Hamiltonian of {Gd2} takes the form:

$$\hat H_{\mathrm{spin}}(\{ {\mathrm{Gd}}_2\} ) = - 2j_{12}\hat S_{{\mathrm{Gd}}_1}\hat S_{{\mathrm{Gd}}_2} - 2\hat s\left( {K_1\hat S_{{\mathrm{Gd}}_1} + K_2\hat S_{{\mathrm{Gd}}_2}} \right) \\ \approx - 2K^{{\mathrm{eff}}}\hat s(\hat S_{{\mathrm{Gd}}_1} + \hat S_{{\mathrm{Gd}}_2}),$$
(4)

where \(\hat S_{{\mathrm{Gd}}_i}\) denotes the Gd spin at site i. Density-functional theory (DFT) calculations predict two very close Ki values for {Gd2}, 181 and 184 cm−1, and much smaller direct Gd–Gd exchange coupling, j12=−1.2 cm−1.38 When j12 is that small, the Ki and j12 parameters cannot be determined from experimental data separately, leading to the approximate Hamiltonian Eq. 4 with an effective coupling Keff (ref. 45). Having used this single-parameter Hamiltonian in the analysis of the measured magnetization and χT curves for {Gd2} we have obtained the Keff value of 160±10 cm−1 (Supplementary Figs 34, 36) which is close to the one of the [Gd3+e–Gd3+] spin-system in Gd2@C79N (Keff = 170±10 cm−1)45 and is in good agreement with the DFT results.

Strong exchange interactions result in the giant-spin magnetic ground state of {Gd2} with S = 15/2, created by two local Gd spins (SGd = 7/2) ferromagnetically coupled via the free electron spin (Se = 1/2). Fine details of this state can be further attested by EPR spectroscopy. At room temperature in toluene solution, {Gd2} shows a single EPR line with a g-factor of 1.987 (Supplementary Figs 37, 38). Freezing the solution at 100 K results in a complex multiline structure in the X-band (9.4 GHz) EPR spectrum (Fig. 3a). The Q-band (34 GHz) spectrum of {Gd2} measured under similar conditions has a simpler but still complex pattern. The low-temperature structure in the EPR spectra is an evidence of zero-field interactions in the large-spin ground state. The spin Hamiltonian of such a system can be written as:

$$\hat H_{\mathrm{spin}} = D\left( {\hat S_z^2-\frac{1}{3}S\left( {S + 1} \right)} \right) + \frac{1}{2}E\left( {\hat S_ + ^2 + \hat S_-^2} \right) + g_{{\mathrm{iso}}}{\mathrm{\mu }}_{\mathrm{B}}{\boldsymbol{B}}\hat S,$$
(5)

where the first two terms describe the second-order zero-field splitting (ZFS) of rhombic symmetry, and the last term represents the Zeeman effect. The X-band and Q-band EPR spectra of {Gd2} in a frozen solution can be well reproduced by the parameters D = 1.00(2) GHz, E = 0.22(4) GHz, and giso = 1.987 (Fig. 3a). The ZFS tensor of {Gd2} is found to be similar to that of the previously reported Gd2@C79N (D = 0.96(6) GHz, E = 0.14(1) GHz, giso = 1.99)44,55, but shows somewhat larger rhombicity, which is in line with the asymmetric geometry of {Gd2} induced by the exohedral CH2Ph group. A schematic description of the Zeeman splitting of the 16 energy levels of the weakly anisotropic S = 15/2 system in {Gd2} together with the transitions accessible in the X- and Q-band EPR spectra are shown in Fig. 3b and Supplementary Fig. 39.

Fig. 3
figure 3

Electron paramagnetic resonance (EPR) spectroscopy of {Gd2}. a X-band and Q-band EPR spectra of frozen {Gd2} solution in toluene near 100 K together with the spectra simulated for spin S = 15/2 with giso = 1.987 and zero field splitting (ZFS) parameters D = 1.00 GHz and E = 0.22 GHz (inhomogeneous broadening is accounted for by ZFS strain StrD = 0.029 GHz and StrE = 0.027 GHz); asterisks mark unidentified signals (presumably of low spin states or organic impurities), the inset shows the spin-density distribution in {Gd2}; b Zeeman splitting for spin S = 15/2 with the above ZFS parameters (magnetic field is parallel to z-axis of the ZFS tensor); also shown are energies of the X-band (9.4 GHz) and Q-band (34 GHz) microwave photons, EPR-active transitions (ovals and small arrows), and the resonance fields corresponding to the g-factor of 1.987 (vertical dotted lines)

Blocking of magnetization in {Ln2}

An essential characteristic of SMM is the blocking of magnetization at a certain temperature (when relaxation of magnetization becomes too slow at the time scale of the measurement). This temperature can be identified by a characteristic divergence of magnetic susceptibilities χFC and χZFC, measured for a field-cooled and a zero-field cooled sample, respectively. The blocking temperature of magnetization, TB, is usually defined as the peak in χZFC. Four {Ln2} compounds exhibit blocking of magnetization above 2 K. The TB value of {Tb2}, 28.9 K (Fig. 4a), is one of the highest among all known SMMs12,18,30,56,57,58. {Dy2} also features a high TB value of 21.9 K38. The χZFC of the mixed Tb-Gd sample shows two peaks at 14.4 K and near 29 K (Fig. 4a). The latter corresponds to {Tb2}, whereas the peak at 14.4 K can be assigned to {TbGd} since {Gd2} does not show blocking of magnetization above 1.8 K. Finally, in {TbY} χFC and χZFC diverge below 5 K, although there is no peak in χZFC. Another universal SMM parameter, the 100-s blocking temperature TB100, is determined from relaxation times of magnetization (see below) to be 18.2 K for {Dy2} and 25.2 K for {Tb2}. TB100 of {Tb2} is surpassed only by a recently reported group of Dy-metallocenium salts with different alkyl groups in cyclopentadienyl rings, which show TB100 values of 53–65 K12,13,14,15.

Fig. 4
figure 4

Magnetic properties of {Ln2} molecules. a Blocking temperature of magnetization in {Tb2}, {TbGd} and {TbY}; dotted lines are measurements of magnetic susceptibility χ during cooling in the field of 0.2 T, solid lines are measurements during heating in the field of 0.2 T of zero-field cooled samples (sweep rate 5 K min−1, arrows indicate direction of the measurement for each curve), vertical red dotted lines denote TB values; b magnetic hysteresis curves for {Tb2}, sweep rate 9.5 mT s−1; c alignment of Ln magnetic moments in {Ln2} according to ab initio calculations: collinear in {Tb2} and {Dy2}, tilted in {Ho2} (the arrows indicate directions of the single-ion quantization axis for each Ho), easy-plane in {Er2}, in the latter the Ln spins are visualized as ellipsoids; d low-temperature magnetization curves for {TbY}, sweep rate 2.9 mT s−1; the inset shows enhancement of the field range, in which magnetic hysteresis is observed

Magnetic properties of {Tb2} and {Dy2}

In accordance with its high TB, {Tb2} shows magnetic hysteresis up to a temperature of 28 K (sweep rate 9.5 mT s−1). The hysteresis is extremely broad (Fig. 4b), with giant coercive fields of 8 T at 10 K and 8.2 T at 5 K. This value is similar to that of the recently reported dinuclear Tb-metallocene with N23− radical bridge30 (which also has a high TB100 of 20 K) and has no further analogs among molecular magnets or bulk magnetic materials. {Dy2} also exhibits magnetic hysteresis up to 21 K with a coercive field of 1.2 T at 1.8 K38.

The relaxation times of magnetization below TB are determined by stretched exponential fitting of magnetization decay curves (Supplementary Figs 40, 41). Between 2 and 15 K in zero field {Tb2} exhibits a temperature-independent relaxation time of (6.5 ± 1) × 104 s, which is an indication of the quantum tunneling of magnetization (QTM, Fig. 5a). In this relaxation regime, the giant combined spin of the [Tb3+e–Tb3+] system flips as a single entity. When QTM is quenched by a finite magnetic field, the relaxation times of {Tb2} are further increased by orders of magnitude. A conservative estimation for the relaxation time in the field of 0.3 T is reaching 6 years at 3 K. Above 20 K the relaxation of magnetization in {Tb2} shows a linear temperature dependence in Arrhenius coordinates. Between 35 and 45 K the relaxation times are determined from AC measurements (Supplementary Fig. 42). The τm values continue the linear regime found in DC measurements, which is well described by the Orbach relaxation mechanism, \(\tau _{\mathrm{m}}^{ - 1} = \tau _0^{ - 1}{\mathrm{exp}}( - U^{{\mathrm{eff}}}/T)\), with Ueff = 799 ± 2 K and τ0 = (1.66 ± 0.14) × 10−12 s. Here Ueff is the effective barrier corresponding to the excited spin state involved in the relaxation, and τ0 is the attempt time. The Ueff value determined for {Tb2} is the largest among all radical-bridged lanthanide molecule magnets described, hitherto.

Fig. 5
figure 5

Relaxation of magnetization in {Ln2} molecules. a magnetic relaxation times of {Tb2}, full dots are zero-field data, open dots are in-field data, red line denotes Orbach processes, solid horizontal line denotes QTM; the inset shows schematically two main relaxation pathways in {Tb2}: QTM and Orbach relaxation via an exchange-excited state, red arrows are Ln spin, blue arrow is a free electron spin; b magnetic relaxation times of {Ho2} at zero field (full dots) and in different fields between 0.1 and 0.4 T, red and purple solid lines are Orbach processes, blue line is a possible Raman contribution (~T10.1); the inset shows magnetic field dependence of relaxation times at 1.8 and 5 K; c magnetic relaxation times of {TbY}, dashed horizontal line is a QTM contribution to zero-field relaxation, magenta line is a low-power process (~T1.70), dark blue line is a combination of both, light blue line is a Raman process (~T 4.64)

The relaxation of magnetization in {Dy2} follows the same trend as in {Tb2}, but with lower temperatures and shorter times38. A zero-field QTM regime with τQTM of (3.3 ± 0.2) × 103 s is found for {Dy2} below 5 K. Above 20 K, relaxation of magnetization in {Dy2} is described by the Orbach mechanism with Ueff = 613 ± 8 K and τ0 = (3.6 ± 1.0) × 10−12 s.

The spin Hamiltonian for anisotropic lanthanides requires single-ion ligand field parameters, which were computed here ab initio. At first, the {Ln2} molecules were optimized at the DFT level (Supplementary Tables 36), and then the LF states were computed at the CASSCF/RASSI-SO level for each Ln center in a model {LnY} system, where Y substitutes one of the lanthanide ions (Supplementary Tables 711). The calculations revealed that Tb3+ and Dy3+ in the {Ln2} molecules have easy-axis anisotropies with a high-spin ground state, which are described by \(| \pm 15/2\rangle\) and \(| \pm 6\rangle\) doublets for Dy and Tb, resp. (complete details of the LF splitting are described in Supplementary Table 8, possible origins of the axial ligand field imposed on lanthanide ions in {Ln2} molecules are discussed in Supplementary Notes 2 and 3, Supplementary Fig. 43, Supplementary Table 12 and in ref. 38). Thus, the ground-state spins of Tb or Dy are of Ising type with their easy axes aligned along the Ln–Ln bond. Therefore the low-energy states of the [Ln3+e–Ln3+] systems can be described by the same effective spin Hamiltonian as in Eq. 4, but with addition of the ligand field terms:

$$\hat H_{{\mathrm{spin}}}(\{ {\mathrm{Ln}}_2\} ) = \hat H_{{\mathrm{LF}}_1} + \hat H_{{\mathrm{LF}}_2} - 2K^{{\mathrm{eff}}}\hat s\left( {\hat J_{{\mathrm{Ln}}_1} + \hat J_{{\mathrm{Ln}}_2}} \right).$$
(6)

Here again, the direct Ln–Ln exchange is neglected. A similar form of the Hamiltonian can be derived from the Lines model59, assuming isotropic coupling for each Kramers Doublet (KD). Although a description of the exchange in Eq. 6 as isotropic is an oversimplification, for the Ising ground state with collinear spins this approximation is essentially valid. Thus, the Hamiltonian Eq. 6 provides a correct description of the ground state properties of the systems with strong easy-axis anisotropy, such as {Tb2} and {Dy2}, but it is expected to be less reliable when interactions with higher KDs are involved. Since the relative energies of the first and second excited doublets, KD2 and KD3, are predicted to be in the range of 300–500 cm−1 for Tb and 200–400 cm−1 for Dy, only the lowest-energy excited states contribute to the magnetic properties in the experimentally relevant temperature range.

The values of Keff for {Tb2} and {Dy2} can be estimated by modeling the susceptibility and the magnetization curves using Hamiltonian Eq. 6 with the addition of the Zeeman term60. The shapes of these curves (Supplementary Figs 4446) correspond to the strong coupling with the Keff values of 45–53 cm−1 for {Tb2} and 30–35 cm−1 for {Dy2}. For such large values, the estimates from χT fitting give rather large confidence limits. Yet, a more precise estimation is possible through modeling of the effective barriers of the Orbach process, Ueff. It can be done assuming that the low-energy LF excited states do not participate in the relaxation of {Tb2} and {Dy2}, while an efficient Orbach process involves the first exchange-excited state corresponding to flipping of one of the lanthanide spins under the Hamiltonian Eq. 6 (Fig. 5a)30,38,49. If only the exchange term of the Hamiltonian Eq. 6 is considered, and the ground state lanthanide spins are of Ising type with Jz = ±J (here J is the total momentum of the lanthanide ion), the energy of the exchange-excited state and hence the relaxation barrier would be Ueff=2JKeff with Keff of 46 cm−1 for {Tb2} and 28 cm−1 for {Dy2} (see ref. 30 for a similar discussion on the radical-bridged dinuclear Tb complex). However, mixing of the LF and exchange excitations changes the energies of the states with predominant exchange excitation (Supplementary Table 13, Supplementary Fig. 47, Supplementary Note 4), resulting in the Keff values of 55 cm−1 and 32 cm−1 for {Tb2} and {Dy2} systems, respectively. For comparison, the largest lanthanide-radical coupling constants in the radical-bridged complexes are −23.1 cm−1 for Tb and −7.2 cm−1 for Dy in [Ln3+–N23−–Ln3+] systems30.

Magnetic properties of {Ho2}

Relaxation times of {Ho2} determined by AC magnetometry between 1.8 and 20 K (Supplementary Figs 4853) show peculiar temperature and field dependences (Fig. 5b). Between 1.8 and 10 K in zero DC field, the compound exhibits a linear log(τm)-vs-T−1 dependence resembling the Orbach mechanism with Ueff = 5.3 ± 0.1 K and τ0 = 2.46 ± 0.06 ms. With the increase of the magnetic field up to 0.5 T, the relaxation decelerates, and the linear dependence is gradually transformed into a curved one (Fig. 5b). If the magnetic field exceeds 0.5 T, the relaxation accelerates again. Deceleration of the relaxation of magnetization with the application of a magnetic field is a typical characteristic of QTM, but the pronounced temperature dependence is not common for a ground-state QTM. However, Zheng and Chilton et al. have recently emphasized that QTM may show temperature dependence due to the temperature-dependent phonon collision rate61, and temperature-dependent QTM was also observed in DySc2N@C8062. Alternatively, such a behavior may correspond to the thermally-assisted QTM, i.e., QTM in an excited state, but the nature of this state is not clear since the energy of 5.3 K is much smaller than the predicted ligand field splitting (see below). The increase of the relaxation rate with the further increase of the magnetic field is an indication of the relaxation via the direct mechanism, which involves phonons of the frequency corresponding to the energy gap between the opposite spin states.

Above 10 K, the relaxation of magnetization in {Ho2} is field-independent and exhibits rapid acceleration with temperature. The log(τm)-vs-T−1 dependence deviates from the linear form and cannot be assigned to a single Orbach process (Fig. 5b). The best fit to the data is obtained by a combination of Orbach and Raman mechanisms (the latter implies a relaxation rate as power function of temperature, \(\tau _{\mathrm{m}}^{ - 1} = AT^n\)) with Ueff = 334 ± 10 K, τ0 = (5.6 ± 2.6) × 10−13 s, n = 10.1 ± 0.3, and A = (1.71 ± 1.3) 10−9 s−1 K−10.1. A similarly good fit is obtained for a combination of two Orbach processes with parameters U1eff = 324 ± 5 K, τ01 = (0.8 ± 0.2) × 10−12 s, and U2eff = 136 ± 5 K, τ02 = (1.1 ± 0.4) × 10−7 s.

Ab initio calculations for Ho3+ in {Ho2} predict a high-spin ground state and smaller energy splitting compared to {Tb2} and {Dy2}. Further, the easy-axis of each ion is tilted from the Ho–Ho axis by 13.4° (Fig. 4c), and the quasi-doublet ground states have strongly mixed mJ character due to higher-order LF terms (Supplementary Tables 9 and 12, Supplementary Note 3): the leading mJ terms are 64% \(| \pm 8\rangle\), 14% \(| \pm 7\rangle\), and 10% \(| \pm 6\rangle\). The lower energy splitting and non-collinearity of the single-ion magnetic moments are consistent with the smaller paramagnetic shifts observed for {Ho2} in the 1H NMR spectra when compared to {Tb2} and {Dy2}.

For non-collinear magnetic moments of Ho ions in {Ho2} the anisotropy of the exchange coupling should have a pronounced effect already in the ground state, and the spin Hamiltonian Eq. 6 cannot capture the whole underlying physics. Indeed, our attempts to reproduce the experimental χT and magnetization curves of {Ho2} using Hamiltonian Eq. 6 gave poorer agreement than for {Tb2} and {Dy2} (Supplementary Figs 5456). The closest match between experiment and simulations is found for Keff = 40 cm−1. With this value of the coupling constant, the lowest-energy exchange-excited state is predicted at 374 K (Supplementary Table 14, Supplementary Fig. 57), which can be compared to the estimated Ueff value of 324–334 K.

The non-collinearity effect can alternatively be introduced via phenomenological Dzyaloshinkii-Moriya interaction term between the lanthanide sites, \(\hat H_{DMI} = {\mathbf{D}}_{12} \cdot (\hat J_{{\mathrm{Ln}}_1} \times \hat J_{{\mathrm{Ln}}_2})\), where D12 is the Dzyaloshinskii–Moriya vector63. Presence of this term leads to staggered magnetization and also allows for the exchange and Kondo spin-fluctuation processes, providing a spin relaxation mechanism for {Ho2}.

Magnetic properties of {Er2}

For {Er2}, dynamic susceptibility measurements showed the out-of-phase χ” response only at low temperature and in the presence of a magnetic field (Supplementary Figs 5860). At 1.8 K, the relaxation time increases with the field from 18 ms at 0.1 T to 47 ms near 0.5 T, but the amplitude of χ” is the highest in the field of 0.25 T. The temperature dependence measured in the field of 0.25 T revealed a decay of the relaxation time from 30 ms at 1.8 K to 10 ms near 3 K; at higher temperatures the signal becomes too weak to be measured reliably. The assignment of the underlying relaxation process to the SMM features of {Er2} is ambiguous because similar relaxation behavior at low temperatures may be also caused by the lattice-based phonon bottleneck. In good agreement with paramagnetic NMR data, our ab initio calculations predict an easy-plane character of the ground state of Er3+ ions in {Er2} as visualized in Fig. 4c (see also Supplementary Table 10 and extended discussion in Supplementary Note 3). In this situation the simple spin Hamiltonian Eq. 6 cannot describe the spin exchange processes well (Supplementary Figs 6163). Moreover, the easy-plane anisotropy implies the presence of strong spin-fluctuation processes, which provides spin relaxation mechanism detrimental to good SMM behavior.

Magnetic properties of {TbY}

To understand, if the symmetry of the [Ln3+e–Ln3+] spin system is essential for the excellent SMM performance of {Tb2}, we studied mixed-metal {TbGd} and {TbY} compounds. The blocking of magnetization near 14 K in {TbGd} shows that coupling the large isotropic spin of Gd to the anisotropic spin of Tb via the delocalized electron spin gives reasonably strong SMM, yet the absence of anisotropy on one of the metal sites leads to the two-fold decrease of the blocking temperature.

Substitution of one Tb ion in {Tb2} by a non-magnetic Y results in a dramatic increase of the relaxation rate. {TbY} shows narrow magnetic hysteresis only below 5 K (Fig. 4d, Supplementary Figs 64, 65). The opening is observed in the field range of 0.1–1.0 T, whereas near zero field the hysteresis loop is closing. AC measurements (Supplementary Figs 6668) also showed that in-field and zero-field magnetization relaxation times of {TbY} are considerably different below 15 K (Fig. 5c), and the difference is reaching a factor of 450 at 2 K (2.9 s at 0.3 T versus 6 ms in zero field). The field dependence of τm measured at 4 K has a sharp maximum at 0.25 T. Such a strong variation of relaxation time with the magnetic field points to a considerable contribution of zero-field QTM at low temperature. However, the zero-field relaxation rate shows temperature dependence down to 1.8 K. The low-T part can be well described by a combination of temperature-independent QTM and a power function of temperature, \(\tau _{\mathrm{m}}^{ - 1} = \tau _{{\mathrm{QTM}}}^{ - 1} + AT^n\), with τQTM = 19.0 ± 0.6 ms, A = 16 ± 1 s−1 Kn, and n = 1.70 ± 0.04. The exponent of 1.7 is close to the expected value for a direct (n = 1) or a bottlenecked direct process (n = 2). However, this temperature-dependent process should be strongly linked to the QTM because it is not observed anymore when the finite field of 0.3 T is applied. Temperature dependence of the in-field relaxation rate as well as high-temperature zero-field relaxation are well described by a power function with parameters A = 2.5 ± 0.5 ms−1 Kn, and n = 4.64 ± 0.08. Fitting of χT and magnetization measurements with the short version of the Hamiltonian Eq. 6, including only single lanthanide ion exchange-coupled to electron spin gives the optimal Keff value of 35 cm−1 (Supplementary Figs 6971), which is considerably smaller than the Tb-electron coupling constant in {Tb2}. These results prove that the coupling of the single lanthanide spin to a delocalized electron spin of the single-electron Tb–Y bond is not sufficient to create a strong SMM and that the presence of two local lanthanide spins in {Ln2}, preferably both of uniaxial anisotropy type, is indeed essential.

Electrochemistry and properties of {Ln2} anions

The Ln–Ln bonding orbital occupied by a single electron is expected to be redox-active. All {Ln2} compounds exhibit reversible electrochemistry (Fig. 6a shows the cyclic voltammogram of {Er2} as a representative example, see Supplementary Figs 7277 for other {Ln2} molecules) with almost identical first oxidation potentials at 0.50−0.52 V and strongly metal-dependent first reduction potentials varying from −0.86 V in {Gd2} to −0.42 V in {Er2} (Fig. 6, potentials are referred versus Fe(Cp)2+/0 redox couple and listed in Supplementary Table 15). The first reduction potentials correlate well with Shannon ionic radii64 of the metals; even better correlation is found between the first reduction potentials and the 4fn5d16s2 → 4fn5d26s1 excitation energies of lanthanide atoms (Supplementary Fig. 78). The metal dependence is a clear indication of the population of the single-electron Ln–Ln bonding orbital by the second electron in the {Ln2} anion as sketched in Fig. 6b. The formal Ln2.5+ oxidation state of the pristine {Ln2} is thus transformed into the Ln2+ state in {Ln2}.

Fig. 6
figure 6

Electron transfer properties of {Ln2} molecules. a Cyclic voltammogram of {Er2} in o-dichlorobenzene solution as a representative example of the {Ln2} series. b Schematic description of the single-electron reduction and oxidation of {Ln2} compounds showing addition of one electron to the Ln–Ln bond and removal of one electron from the fullerene cage. c The first oxidation (red dots) and reduction (blue dots) potentials of {Ln2} in o-dichlorobenzene/TBABF4 solution as a function of ionic radius of Ln (for {TbY}, the average radius of Tb3+ and Y3+ is used; lines are shown to guide the eye); d 1H NMR spectra of {Tb2}, {Ho2} and {Er2} anions in d4-o-dichlorobenzene (colored lines) in comparison to the spectra of neutral compounds (light gray lines). e Schematic description of the spin-valve effect of the {Ln2} molecule: in a certain bias range limiting the current to the metal-based LUMO, only the electrons with their spin antiparallel to the spin of the molecule can pass through

Transformation of the single-electron Ln–Ln bond into the two-electron covalent bond upon the first reduction should decrease the exchange coupling between the lanthanide spins. The system of two weakly-coupled spins in {Ln2} would have a much smaller magnetic anisotropy and hence a smaller paramagnetic shift of the 1H NMR resonances than the strongly coupled dimer in {Ln2}. Indeed, 1H NMR spectra of the anions {Tb2} and {Ho2}, produced by a single-electron reduction of pristine {Ln2} compounds with cobaltocene, show a considerable decrease of the 1H paramagnetic shifts (Fig. 6d). Surprisingly, in {Er2} the dipolar shifts are increased when compared to the pristine {Er2}. The latter indicates that the anisotropic exchange interactions in {Er2} are enhancing the axial (z) component of the magnetic susceptibility and hence decrease the total magnetic anisotropy in comparison to the weakly-coupled spin system in {Er2}.

Discussion

To summarize, we showed that just like delocalization of electrons is the essence of covalent chemical bonding, delocalization of the unpaired electron spin in the [Ln3+e–Ln3+] system between two lanthanide sites glues their magnetic moments together. The strong ferromagnetic coupling emerging from these interactions is responsible for the high-spin magnetic ground state in certain {Ln2} molecules. But the strong exchange coupling alone is not enough to make a good SMM. Single-ion anisotropy and collinearity of lanthanide spins play a crucial role as well, and whereas {Tb2} and {Dy2} with high-spin easy-axis single-ion ground states and collinear moments make the best SMMs in the series, {Ho2} with mixed LF states and tilted magnetic moments is only a modest SMM, and {Er2} with easy-plane single-ion anisotropy is hardly an SMM at all. The different magnetic anisotropy of these chemically very similar molecules can be understood from the element-specific shape of the Ln-4f charge density, interacting with the molecular charge distribution and being spin-orbit coupled to the 4f magnetic moments (Supplementary Note 3). Furthermore, the coupling of a single lanthanide spin to the delocalized electron spin is also not sufficient as illustrated by the soft SMM behavior of the mixed-metal {TbY} system. Thus, homonuclear lanthanide dimers with collinear magnetic moments and strongly-axial single-ion magnetic anisotropy give the best SMMs.

Importantly, although the Ln2 dimer is protected by the fullerene, it is not completely isolated from the environment. The carbon cage remains transparent for electrons65,66, and {Ln2} compounds exhibit lanthanide-based redox-activity. In the first reduction step, the Ln–Ln bonding orbital is populated by a second electron, thus allowing to change the valence state from Ln+2.5 to Ln+2. Simultaneously, the exchange interactions between Ln spins are reduced in the anionic state. Importantly, the valence electrons are strongly coupled to the lanthanide spins. Thus, magnetic properties of {Ln2} molecules can be switched by an electron transfer, which forms a background for their possible application in spin-polarized molecular transport, as redox magnetic switches or as an electron spin detector using magnetoresistance (Fig. 6e).

Methods

Synthesis

{Ln2} compounds were produced by the Krätschmer-Huffman method followed by functionalization with a benzyl group using the method developed in ref. 38 (Supplementary Fig. 1) The graphite rods (length 100 mm, diameter 6 mm) are packed with metal oxides mixed with graphite (molar ratio of Ln:C = 1:15) and evaporated in an electric arc with a current of 100 A in 180 mbar helium atmosphere. The fullerene-containing soot is extracted under nitrogen for 20 h by boiling dimethylformamide (DMF), and DMF solution of EMFs is then reacted with excess of benzyl bromide BrCH2Ph for 20 h at elevated temperature under nitrogen protection. Afterwards the solvent is evaporated, and the residue is washed with methanol. The rest is dissolved in toluene and further separated by high performance liquid chromatography with Buckyprep, and Buckyprep-D columns (Nacalai Tesque, Japan) as shown in Supplementary Figs 219. 0.5–2 mg of pure {Ln2} compounds could be isolated in this work. The yield of {Ln2} depends on the metal size: the yields of {Er2}, {Ho2}, and {Dy2} are similar, the yield of {Tb2} is ca twice lower, and the yield of {Gd2} is the lowest in the series.

Spectroscopic and electrochemical measurements

Matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) mass-spectra were measured with a Bruker autoflex mass-spectrometer with 1,1,4,4-tetraphenyl-1,3-butadiene as a matrix. EPR spectra of {Gd2} solution in toluene were measured using cw-EPR spectrometer EMX Plus (Bruker), working in X-band and Q-band regions. The EPR spectra were fitted using Easyspin, a MATLAB toolbox67. UV-vis-NIR absorption spectra were measured in toluene solution at room temperature with a Shimadzu 3100 spectrophotometer. Raman spectra were recorded at 78 K on a T 64000 triple spectrometer (Jobin Yvon) using 656 nm excitation wavelength of a tunable dye laser Matisse 2 (Sirah Lasertechnik). IR spectra were measured at room temperature with Vertex 80 FTIR spectrometer (Bruker) equipped with a Hyperion microscope. For Raman and IR measurements, the {Ln2} samples were drop-casted from toluene solution onto single-crystal KBr disks. NMR spectra were measured with an Avance 500 spectrometer (Bruker). Voltammetric experiments were performed with a sweep rate of 100 mV s−1 in o-dichlorobenzene solution with TBABF4 electrolyte salt in an oxygen-free glove box using potentiostat-galvanostat PARSTAT 4000 A. A three-electrode system with a platinum working and a counter electrode and a silver wire pseudo-reference electrode was used, potentials were calibrated by adding ferrocene as an internal standard in the end of each measurement.

Single-crystal X-ray diffractometry

Crystal growth of Dy2@C80-CH2Ph۰0.67(CH3Ph) was accomplished by layering hexane over a solution of {Dy2} in toluene. Slow diffusion of two solutions resulted in formation of small black crystals (30 × 30 × 10 μm3). X-ray diffraction data have been collected at 100, 130, 160, 190, 220, 250, 270, and 290 K on the BL14.3 beamline operated by the Joint Berlin MX Laboratory at the BESSY II electron storage ring (Berlin-Adlershof, Germany)68 using a MAR225 CCD detector, λ = 0.89429 Å. Processing diffraction data was done with XDSAPP2.0 suite69. The structure was solved by direct methods and refined using all data (based on F2) by SHELX 201670. Hydrogen atoms were located in a difference map, added geometrically, and refined with a riding model. Crystal data and data collection parameters are summarized in Supplementary Table 1

Magnetometry

DC magnetic measurements were performed using a Quantum Design VSM MPMS3 magnetometer. The samples were drop-casted from CS2 solution into a standard powder sample holder. Long magnetization relaxation times of {Tb2} were determined from the measurement of magnetization decay using a dc-SQUID. After the sample was magnetized to the saturation at 7 Tesla, the field was swept fast to zero or 0.3 T, and then the decay of magnetization was followed over hours and fitted with stretched exponential (Supplementary Figs 40, 41). Measurements of magnetic hysteresis curves of {Tb2} were accomplished with a PPMS system equipped with a 14 T magnet. AC-magnetometry measurements were performed using Quantum Design MPMS XL magnetometer, Quantum Design VSM MPMS3 magnetometer, and PPMS system for the high-frequency range (0.5–10 kHz). See Supplementary Methods for further details.

Calculations

VASP code, version 5.0, was employed to optimize the molecular structures at the PBE-D level using PAW pseudopotentials71. Ab initio energies and wave functions of CF multiplets for the {LnY} molecules have been calculated at the CASSCF/SO-RASSI level of theory using the quantum chemistry package MOLCAS 8.072. The single ion LF-parameters were calculated based on ab initio data with the use of SINGLE_ANISO module73. Modelling of the magnetic properties was accomplished with the PHI program60 and included powder-averaging. See Supplementary Methods for further details.