A Trinuclear Gadolinium Cluster with a Three-Center One-Electron Bond and an S = 11 Ground State

The recent discovery of metal–metal bonding and valence delocalization in the dilanthanide complexes (CpiPr5)2Ln2I3 (CpiPr5 = pentaisopropylcyclopentadienyl; Ln = Y, Gd, Tb, Dy) opened up the prospect of harnessing the 4fn5dz21 electron configurations of non-traditional divalent lanthanide ions to access molecules with novel bonding motifs and magnetism. Here, we report the trinuclear mixed-valence clusters (CpiPr5)3Ln3H3I2 (1-Ln, Ln = Y, Gd), which were synthesized via potassium graphite reduction of the trivalent clusters (CpiPr5)3Ln3H3I3. Structural, computational, and spectroscopic analyses support valence delocalization in 1-Ln resulting from a three-center, one-electron σ bond formed from the 4dz2 and 5dz2 orbitals on Y and Gd, respectively. Dc magnetic susceptibility data obtained for 1-Gd reveal that valence delocalization engenders strong parallel alignment of the σ-bonding electron and the 4f electrons of each gadolinium center to afford a high-spin ground state of S = 11. Notably, this represents the first clear instance of metal–metal bonding in a molecular trilanthanide complex, and the large spin–spin exchange constant of J = 168(1) cm–1 determined for 1-Gd is only the second largest coupling constant characterized to date for a molecular lanthanide compound.


Synthesis of (Cp iPr5 ) 3 Gd 3 H 3 I 2 (1-Gd).
Under argon, (Cp iPr5 )2Gd2I4 (0.650 g / 0.473 mmol) and benzene (50 mL) were combined in a 100 mL Schlenk flask to give a pale-yellow solution. Solid LiCH2SiMe3 (0.089 g / 0.95 mmol) was added and the reaction mixture was stirred at room temperature for 24 h, during which time the formation of a yellow/amber suspension with colorless precipitate was noted. Solvent was removed under vacuum to leave a yellow/amber solid residue, which was slurried by stirring with n-hexane (30 mL) for 30 min, then filtered through Celite. The filter pad was extracted with additional n-hexane (2 × 5 mL) to give a yellow/amber filtrate. This yellow-amber solution containing crude (Cp iPr5 )2Gd2(CH2SiMe3)2I2 was added to a Fisher-Porter tube apparatus and frozen with liquid N2 while applying dynamic vacuum. While the solution was still frozen, H2 gas was dosed into the tube at a pressure of 80 psi for approximately 30 s, the apparatus was then sealed, and the reaction mixture was allowed to thaw and stirred at room temperature for 24 h to give a yellow/amber solution Solvent was removed under vacuum to leave a yellow/amber sticky residue containing crude (Cp iPr5 )3Gd3H3I3, which was taken up in n-hexane (50 mL) and filtered through Celite; the pad washed with additional n-hexane (2 × 10 mL) to give a yellow/amber filtrate. This filtrate was transferred to a 100 mL Schlenk flask, KC8 (0.213 g / 1.58 mmol) was added, and the mixture stirred at room temperature for 9 d. The reaction mixture was filtered through Celite to give a dark blue filtrate, and the residue remaining on the filter pad extracted with additional n-hexane (200 mL) and filtered through Celite. The filtrates were combined, concentrated to 60 mL, heated to boiling to re-dissolve a blue precipitate, left to cool to room temperature overnight, then transferred to the freezer (−35 °C). Small dark green-blue, thin rectangular plank-like crystals of 1-Gd were isolated in multiple crops, washed with a small amount of cold (−35 °C) pentane and dried under vacuum (0.120 g / 0.077 mmol / 24% based on (Cp iPr5 )2Gd2I4). A separate batch of 1-Gd crystals was also grown from a concentrated methylcyclohexane solution in the freezer (

X-ray Crystallography Data Collection and Refinement Details
Structure solutions were determined by SHELXT using direct methods and these solutions were refined via least-square refinement against F 2 by SHELXL, as implemented in OLEX2 crystallographic software. 4,5 For all structures, all non-hydrogen atoms were refined anisotropically. All hydrogen atoms except hydride ligands were placed on geometrically calculated positions using the riding model and refined isotropically.
(Cp iPr5 ) 2 Y 2 (CH 2 SiMe 3 ) 2 I 2 . The compound crystallized in the space group P1 with one molecule in the asymmetric unit. The Cp iPr5 ligands showed positional disorder over two sites, and the occupancies of these two components were refined while constraining the sum to unity, yielding a ratio of 0.501(5):0.499 (5) and 0.727 (10):0.273 (10), respectively. The application of the standard combination of SIMU, DELU restraints and EADP constraints to the displacement parameters of the cyclopentadienyl carbons-and EADP constraints to the displacement parameters of the analogous carbon atoms bound to the Y atom in the CH2SiMe3-groups gave free R factors of R1 = 3.38% and wR2 = 7.85%. In checkcif, one alert B (PLAT230_ALERT_2_B) was generated from the metal bound carbon atoms in the structure. No A alerts were found.
(Cp iPr5 ) 2 Gd 2 (CH 2 SiMe 3 ) 2 I 2 . The compound crystallized in the space group P1 with half of a molecule in the asymmetric unit and an inversion center located in the middle of two Gd atoms within the molecule. The Cp iPr5 ligand, the methylene carbon atom in the CH2SiMe3 group, and the iodide ligand showed positional disorder over two sites, and the occupancies of these two components were refined while constraining the sum to unity, yielding ratios of 0.576 (10):0.424 (10), 0.52(3):0.48 (3), and 0.9446 (17):0.0554 (17), respectively. The application of the standard combination of SIMU, DELU restraints to the displacement parameters of the cyclopentadienyl carbons, and EADP constraints to the displacement parameters of the disordered carbon atom in the CH2SiMe3 group and the disordered iodine atom gave free R factors of R1 = 3.10% and wR2 = 7.07%. In checkcif, two level B alerts (PLAT213_ALERT_2_B and PLAT230_ALERT_2_B) were generated from the metal bound carbon atoms in the structure.

(Cp iPr5 ) 3 Y 3 H 3 I 3 .
The compound crystallized in the space group P21/n with one molecule in the asymmetric unit. Solvent mask, as implemented in Olex2 (analogous to SQUEEZE), was used to remove the electron density from the highly disordered solvent molecules in the lattice (V = 1754 Å 3 , 300 e − ). One Cp iPr5 ligand and the three iodide ligands showed positional disorder over two sites, and the occupancies of these two components were refined while constraining the sum to unity, yielding a ratio of 0.673(9):0.327(9) and 0.9608(7):0.0392(7), respectively. The hydride ligands near the Y atoms were located manually using distance restraints (DFIX or SADI). The application of EADP constraints to the displacement parameters of the cyclopentadienyl carbons and iodide ligands gave free R factors of R1 = 4.84% and wR2 = 13.44%. In checkcif, four level B alerts (PLAT971_ALERT_2_B) were generated from positive residual density near the iodide ligands. This residual density could not be modelled as any chemically reasonable species, and it likely arises due to strongly absorbing Y or I or due to problems with the absorption correction.

(Cp iPr5 ) 3 Gd 3 H 3 I 3 .
The compound crystallized in the space group P21/n with one molecule in the asymmetric unit. Solvent mask, as implemented in Olex2 (analogous to SQUEEZE), was used to remove the electron density from the highly disordered solvent molecules in the lattice (V = 1736 Å 3 , 1018 e − ). One Cp iPr5 ligand and iodide ligands showed positional disorder over two sites, and the occupancies of these two components were refined while constraining the sum to unity, yielding a ratio of 0.730(5):0.270(5) and 0.9678(6):0.0322(6), respectively. The hydride ligands near the Y atoms were located manually using distance restraints (DFIX or SADI). The application of the standard combination of SIMU, DELU restraints and EADP constraints to the displacement parameters of the cyclopentadienyl carbons and iodide ligands gave free R factors of R1 = 3.25% and wR2 = 9.05%. In checkcif, one level B alert (PLAT971_ALERT_2_B) were generated from positive residual density near the iodide ligands. This residual density could not be modelled as any chemically reasonable species and likely arises due to strongly absorbing Gd or I or due to problems with the absorption correction.
(Cp iPr5 ) 3 Y 3 H 3 I 2 (1-Y). The compound crystallized in the space group Pbca with one molecule in the asymmetric unit. The structural refinement gives free R factors of R1 = 5.31% and wR2 = 13.50%. In checkcif, one level B alert (PLAT971_ALERT_2_B) was generated from positive residual density near the iodide ligands. This residual density could not be modelled as any chemically reasonable species, and it likely arises due to the strongly absorbing Y or I ions or due to problems with the absorption correction.
(Cp iPr5 ) 3 Y 3 D 3 I 2 (1-Y D ). The compound crystallized in the space group Pbca with one molecule in the asymmetric unit. The structural refinement gives free R factors of R1 = 4.73% and wR2 = 12.44%. In checkcif, three level B alerts (PLAT971_ALERT_2_B and PLAT972_ALERT_2_B) were generated from positive residual density near the iodide ligands. This residual density could not be modelled as any chemically reasonable species and likely arises due to the strongly absorbing Y or I ions or due to problems with the absorption correction.

(Cp iPr5 ) 3 Gd 3 H 3 I 2 (1-Gd).
When crystallized from n-hexane, the compound crystallized in the space group Pbca with one molecule in the asymmetric unit. The structural refinement gave free R factors of R1 = 8.10% and wR2 = 24.09%. In checkcif, three level A alerts (PLAT972_ALERT_2_A) and four level B alerts (PLAT971_ALERT_2_B) were generated from residual density near the iodide ligands and Gd ions. This residual density could not be modelled as any chemically reasonable species, and it likely arises due to the strongly absorbing Gd or iodide ions or due to problems with the absorption correction. In addition, one level B alert (PLAT029_ALERT_3_B) was found due to low data completeness.
In order to resolve these issues, we analyzed crystals of 1-Gd prepared using a different crystallization approach, as described in Section 1. In the second structure obtained from these crystals, the compound crystallized in the space group P1 with one molecule in the asymmetric unit. Solvent mask, as implemented in Olex2 (analogous to SQUEEZE), was used to remove the electron density from the highly disordered solvent molecules in the lattice (V = 232 Å 3 , 53 e − ). The structural refinement gave free R factors of R1 = 4.74% and wR2 = 13.27%. In checkcif, one level B alert (PLAT971_ALERT_2_B) was generated from residual density near the iodide ligands, likely due to the strongly absorbing iodide ions or due to problems with the absorption correction. For consistency with 1-Y, we used the parameters from the first structure in the manuscript.

[(Cp iPr5 ) 3 Gd 3 H 3 I 2 ]
[B(C 6 F 5 ) 4 ] (1-Gd + ). The compound crystallized in the space group Cc with one [(Cp iPr5 )3Gd3H3I2] + cation and one [B(C6F5)4] − anion in the asymmetric unit. Solvent mask, as implemented in Olex2 (analogous to SQUEEZE), was used to remove the electron density from the highly disordered solvent molecules in the lattice (V = 3102 Å 3 , 1425 e − ). Two Cp iPr5 ligands showed positional disorder over two sites, and the occupancies of these two components were refined while constraining the sum to unity, yielding a ratio of 0.759 (14):0.241 (14) and 0.579 (14):0.421 (14), respectively. The hydride ligands near the Gd atoms were located manually using several restraints (DFIX, SADI or FLAT). The application of the standard combination of SIMU, DELU restraints and EADP constraints to the displacement parameters of the cyclopentadienyl carbons and iodide ligands gave free R factors of R1 = 4.34% and wR2 = 11.45%. In checkcif, one level B alert (PLAT342_ALERT_3_B) was generated from the severe disorder within the Cp iPr5 ligands and the tetrakis(pentafluorophenyl)borate counter anion.      Figure S12. Solid-state structure of (Cp iPr5 )2Y2(CH2SiMe3)2I2 with thermal ellipsoids at the 50% probability level. Yellow, purple, turquoise, and grey ellipsoids represent Y, I, Si, and C atoms, respectively. Hydrogen atoms are omitted for clarity. Figure S13. Solid-state structure of (Cp iPr5 )2Gd2(CH2SiMe3)2I2 with thermal ellipsoids at the 50% probability level. Green, purple, turquoise, and grey ellipsoids represent Gd, I, Si, and C atoms, respectively. Hydrogen atoms are omitted for clarity. Figure S14. Solid-state structure of (Cp iPr5 )3Y3H3I3 with thermal ellipsoids at the 50% probability level. Yellow, purple, and grey ellipsoids represent Y, I, and C atoms, respectively. Hydride ligands bound to Y are indicated with white spheres. Hydrogen atoms and the positional disorder of the Cp iPr5 ligand are omitted for clarity. Figure S15. Solid-state structure of (Cp iPr5 )3Gd3H3I3 with thermal ellipsoids at the 50% probability level. Green, purple, and grey ellipsoids represent Gd, I, and C atoms, respectively. Hydride ligands bound to Gd are indicated with white spheres. Hydrogen atoms and positional disorder of the Cp iPr5 ligand are omitted for clarity Figure S16. Solid-state structure of 1-Y with thermal ellipsoids at the 50% probability level. Yellow, purple, and grey ellipsoids represent Y, I, and C atoms, respectively. Hydride ligands bound to Y are indicated with white spheres. Hydrogen atoms are omitted for clarity Figure S17. Solid-state structure of 1-Y D with thermal ellipsoids at the 50% probability level. Yellow, purple, and grey ellipsoids represent Y, I, and C atoms, respectively. Deuteride ligands bound to Y are indicated with white spheres. Hydrogen atoms are omitted for clarity. Figure S18. Solid-state structure of 1-Gd with thermal ellipsoids at the 50% probability level. Green, purple, and grey ellipsoids represent Gd, I, and C atoms, respectively. Hydride ligands bound to Gd are indicated with white spheres. Hydrogen atoms are omitted for clarity. Figure S19. Solid-state structure of 1-Gd + with thermal ellipsoids at the 50% probability level. Green, purple, and grey ellipsoids represent Gd, I, and C atoms, respectively. Hydride ligands bound to Gd are indicated with white spheres. Hydrogen atoms, positional disorder of the Cp iPr5 ligand, and the B(C6F5)4 anion are omitted for clarity.  (11) 3.1435 (7) Ln1 I2 3.1995 (12) 3.1624 (12) 3.1720 (7) Ln2 I1 3.2099 (12) 3.1714 (11) 3.1725 (7) Ln2 I2 3.1500 (12) 3.1622 (12) 3.1323 (7) Ln3 I1 3.1752 (13) 3.1843 (11) 3.1241 (7) Ln3 I2 3.1917 (13) 3.1745 (10) 3.1825 (7) Ln−I (avg) 3.1803 (5)

UV-Vis-NIR Spectra for 1-Y, 1-Gd, and 1-Gd +
UV-Vis-NIR spectra were collected on solutions of 1-Ln (Ln = Y, Gd) in n-hexane and 1-Gd + in 1,2-difluorobenzene. Plots of concentration versus absorbance were used to extract extinction coefficients for each of the features observed in the spectra. Diffuse reflectance spectra were collected for polycrystalline samples of 1-Ln (Ln = Y, Gd) and 1-Gd + , which were diluted in BaSO4 and ground with a mortar and pestle to produce a homogenous powder.

Magnetization and Dc Magnetic Susceptibility Data for 1-Gd and 1-Gd +
Using a Quantum Design MPMS2 SQUID magnetometer, dc magnetic susceptibility data were collected at temperatures ranging from 2 to 300 K under applied fields of 1, 5, and 10 kOe for 1-Gd and 1-Gd + . Magnetic samples were prepared by adding polycrystalline powder (19.0 mg for 1-Gd, 8.2 mg for 1-Gd + ) to a 5 mm i.d./7 mm o.d. quartz tube with a raised quartz platform. A layer of solid eicosane was added on top of each sample (19.0 mg for 1-Gd, 20.2 mg for 1-Gd + ). The tubes were then fitted with Teflon sealable adapters, evacuated using a glovebox vacuum pump, and removed from the glovebox. The portion of the tube containing the sample was cooled quickly in liquid nitrogen, and the tube was flame sealed with an O2/H2 flame. After flame-sealing, the eicosane was melted in a 45 °C water bath in order to provide good thermal contact and to prevent crystallite torqueing. All data were corrected for diamagnetic contributions from the core diamagnetism of the sample and eicosane, estimated using Pascal's constants: χdia = −0.00082865 emu/mol for 1-Gd, −0.00111744 emu/mol for 1-Gd + , and 0.00024306 emu/mol for eicosane.
The dc susceptibility data obtained for 1-Gd and 1-Gd + were fit with the program PHI. 6 Fits to the dc susceptibility data for 1-Gd were performed using the Hamiltonian: where ̂G d1 , ̂G d2 , and ̂G d3 are the spin operators for each metal site in the trinuclear complex and ̂ is the spin operator for the σ-bonding electron, and J is the spin-spin coupling constant.
We have not included Gd-Gd exchange coupling, as there is no sensitivity in the experimental data to fit this value. The only datum we have is the slope of χMT at high temperatures, and here Gd-Gd and 4f-σ parameters are correlated. Thus, we have focused on estimating the strength of the most significant term.
For 1-Gd + , the data were fit using the Hamiltonian where ̂G d1 , ̂G d2 , and ̂G d3 are the spin operators for each metal site and J represents the magnetic exchange coupling between Gd 3+ ions. Figure S34. Zero-field-cooled dc magnetic susceptibility data for 1-Gd under applied fields of 1, 5, and 10 kOe. Black lines represent fits to the data, as described above. Figure S35. Magnetic susceptibility data for 1-Gd in the low temperature region (2 to 50 K). The solid lines correspond to simulations of the data assuming population of only an S = 11 state with g = 1.97, zero field splitting (D = 4.2(2) cm −1 ), and Zeeman splitting using PHI. 6 Figure S36. Zero-field-cooled dc magnetic susceptibility data for 1-Gd + under applied fields of 1, 5, and 10 kOe. Black lines represent fits to the data, as described above. Table S7. Parameters used to fit the dc susceptibility data for 1-Gd and 1-Gd + .

Sample Preparations and Instrument Details
EPR samples were prepared from solutions of (Cp iPr5 )3Y3H3I2 (1-Y) or (Cp iPr5 )3Y3D3I2 (1-Y D ) dissolved in methylcyclohexane solution (ca. 1 mM). The sample used for X-band and Q-band measurements was prepared as follows: the solution was loaded into a quartz tube (o.d. 2.40 mm, i.d. 2.00 mm) inside an Ar-filled glovebox, capped with a rubber septum, and removed from the glove box where the sample was immersed in liquid nitrogen. The sample was then flame sealed under static reduced pressure. The D-band EPR sample was prepared by loading the solution sample into a quartz tube (o.d. 0.600 mm, i.d. 0.500 mm) inside of an Ar-filled glovebox. The sample tube was put inside a 2 mL cryovial, which was capped, then was taken out of the glovebox and immediately submerged in liquid nitrogen. The cryovial was kept submerged in liquid N2 for ~2 h while being transported for measurement (note the cryovial is not airtight). The sample was then removed from the cryovial and directly exposed to air while being loaded into the spectrometer, which had pre-cooled at 90 K under a flow of helium. The sample was kept frozen throughout the whole process.
Note, EPR spectra were not collected for 1-Gd. The first reason for this is that the compound has a non-Kramers ground state and is therefore not necessarily EPR active. For this reason, 1-Y offers a simpler way of studying electron delocalization using EPR spectroscopy. Second, even though EPR signal might be detected, the distribution of zero-field splitting interaction for Gd ions tends to broaden EPR spectra (at high frequencies), obscuring important features, or gives rise to multiple splittings (at X-band frequency), which can be very complex. Lastly, given that the majority (70%) of the isotopes of Gd nuclei are non-magnetic, the use of hyperfine spectroscopies (HYSCORE, ENDOR, EDNMR spectroscopies) to determine the degree of electron delocalization is much more challenging and the results are not as useful as those obtained for 1-Y.
EPR spectra were recorded at the CalEPR center in the Department of Chemistry, University of California, Davis. Continuous wave (CW) X-Band (9.39 GHz) spectra were collected using a Bruker Biospin EleXsys E500 spectrometer (Billerica, MA) equipped with a super high Q resonator (ER4122SHQE). Cryogenic temperatures were controlled and maintained through the use of an ESR900 liquid helium cryostat in conjunction with a temperature controller (Oxford Instruments ITC503) and a gas flow controller. Other spectrometer settings are given in the corresponding figure captions.
Echo-detected field sweep spectra at D-band (129.996 GHz) were collected using a home-built 130 GHz EPR spectrometer equipped with an Oxford-CF935 liquid helium cryostat as described previously. 7 Magnetic field-swept echo-detected EPR spectrum was acquired using the Hahn echo pulse sequence: π/2 − τ − π − echo, with π/2 pulse duration = 35 ns and τ = 250 ns at various S35 magnetic field values. ELDOR-detected NMR (EDNMR) spectra were collected using the pulse sequence The length of the high turning angle pulse (HTA) was set to 14 µs, and the length of microwave pulse /2 was set to 35 ns. The interpulse delays T and τ were set to 6 µs and 250 ns, respectively. A Lorentzian function was used to model the central blind spot and subtracted from the raw EDNMR spectra as previously described. 8 EDNMR spectra were simulated in the same way as ENDOR spectra.
Simulations of all EPR spectra were generated using the Easyspin 6.0.0 toolbox 9 in the Matlab R2021b software suite (Mathworks Inc., Natick, MA).  , and total simulation, respectively. Simulation of Q-band data used 2% impurity, while that of D-band data used 15% impurity. The impurity in the D-band spectrum is attributed to the fact that sample was transported to the EPR facility within an open tube inside of a sealed cryovial, which is not air-tight, and the sample was directly exposed to air while being loaded into the EPR spectrometer, as described in Section 7.1 above. However, the D-band frequency is high enough to separate most of the 1-Y signal (g = [1.9957, 1.9957, 1.965]) from that for the impurity (g = 1.9657), and therefore the impurity is not expected to greatly affect the interpretation of this data. Additionally, the EDNMR spectra collected on the D-band sample (see Figure S41) are impacted by the impurity only at the highest field (4.725 T), and thus the conclusions drawn from these spectra are not meaningfully impacted by the impurity.

EDNMR Spectroscopy
Electron-electron double resonance (ELDOR) spectroscopy is a pulse EPR technique that utilizes two microwave pulses. ELDOR-detected nuclear magnetic resonance (EDNMR) spectroscopy is a technique to study hyperfine interactions of a paramagnetic system. It involves application of a long microwave pulse-termed a high-turning-angle (HTA) pulse-to excite the system from the ground state to (forbidden) excited states, effectively depopulating the ground state. The second pulse is the EPR probe pulse, which typically is the Hanh echo sequence (π/2 -τ -π -τecho) that measures the EPR transition of the system, after it has been altered by the HTA pulse. Figure 43a below depicts a standard Hanh echo sequence used to measure EPR signal, which appears as an echo signal. Figure 43b depicts a typical EDNMR sequence, where a HTA pulse is applied before the Hanh echo sequence. The HTA pulse alters the ground state population of the system and therefore the echo intensity of the EPR signal. By varying the frequency of the HTA pulse (ν), the hyperfine energy levels can be probed.
Consider a simplified example of 1-Y with an unpaired electron (S = 1 /2) coupled to a 127 I nucleus (I = 5 /2). The energy levels of the system can be described by a Hamiltonian: where the first, second, and third terms are electronic Zeeman, electron-nuclear hyperfine, and nuclear Zeeman interactions, respectively (the nuclear quadrupole interaction is ignored for the sake of simplicity). 10 Here, S and I are electron and nuclear spin operators, A is the hyperfine tensor, g is electronic g tensor, gN is nuclear g factor, µB is Bohr magneton, µN is nuclear magneton, and B is magnetic field. Figure 43c displays an energy level diagram of the system upon introduction of the three interactions in respective order. The electronic Zeeman energy is generally much larger than hyperfine and nuclear Zeeman interactions. On the other hand, hyperfine and nuclear Zeeman interactions can be of the same order of magnitude, depending on magnetic field and hyperfine interaction strength. Figure 43c depicts the case of strong coupling (the hyperfine interaction is larger than the nuclear Zeeman interaction).
The black double-headed arrow represents a standard EPR transition (selection rule ΔMS = ±1, ΔMI = 0). Red arrows represent forbidden transitions facilitated by the HTA pulse (ΔMS = ±1, ΔMI = ±1, ±2, …). Note that only a few selected transitions are shown for the sake of simplicity. For full details on EDNMR, see ref. 11). By sweeping the frequency of the HTA pulse (ν), the energy levels of the system may be mapped. To obtain information on nuclear magnetic resonance (NMR) transition (ΔMI = ±1, ΔMS = 0, labeled as blue double-headed arrow in Figure 43c), one subtracts the frequency of EPR transition (ν0) from HTA pulse frequency (ν). The EDNMR spectrum reports signal intensity as a function of |ν−ν0|. The signals appear as negative peaks, as they correspond to the loss of EPR signal, but are subsequently converted to positive peaks by inverting the data. Figure S44 (left) depicts the resulting EDNMR spectra of 1-Y measured at six different field positions along the EPR signal envelope shown on the right. The signal at ν -ν0 = 0 corresponds to a normal EPR transition and is subtracted out through background subtraction (using a Lorentzian function centered at ν -ν0 = 0), generating a blind spot where any data near this position may be concealed. Multiple quantum transitions were also observed and simulated ( Figure S44, orange, yellow, green, and purple lines). These transitions originate from the ΔMI = ±2, ±3, ±4, ±5 transitions (Figure 43c, orange, yellow, green, and purple arrows). Simulations of EDNMR spectra were performed in a similar fashion to that described by Cox et al. 12 The 127 I nucleus was simulated using the salt function of the Easyspin 6.0.0 toolbox. 9 Multiple-quantum transitions were simulated by scaling the frequency relative to the single quantum transitions. The linewidth and intensity of the transitions were manually adjusted to fit the spectra. Black rectangles indicate microwave pulses at the EPR frequency (ν0). Red rectangle represents the HTA pulse, which was applied at a variable frequency ν. (c) Energy level diagram of an S = 1 /2 spin system subjected to electronic Zeeman, electron-nuclear hyperfine, and nuclear Zeeman interactions, respectively. The EPR, EDNMR, and NMR transitions are labeled with black, red, and blue double-headed arrows, respectively. Each state is labeled by the two quantum numbers (MS, MI) to the right. Four additional energy diagrams illustrate the multiple-quantum NMR transitions ΔMI = ±2, ±3, ±4, ±5 (orange, yellow, green, and purple arrows, respectively).

Analysis of Hyperfine Coupling Parameters
Hyperfine interactions may be decomposed into isotropic (Aiso) and anisotropic (T) parts. . The component T is composed of local contribution (Tlocal), which originates from spin density on p, d, or f orbitals, and non-local contribution (Tnon-local), which originates from through-space dipolar interactions:

= local + non−local
Assuming that electron density is equally divided among the three Y ions, we approximate Tnonlocal contributions of ~0.4 and ~0.04 MHz for dipolar interaction between a point charge on three Y nuclei, each with 1 /3 of an electron charge, and the nuclear spin of 127 I and 89 Y, respectively, using the equation: where μ0 is the magnetic permeability of free space and ⃑ ⃑ is a vector connecting the point charge on each Y nucleus with the nucleus of interest. 13 These values are negligible compared to experimental T ( 127 I) and T ( 89 Y) and thus are omitted in the subsequent calculations. An estimate Tnon-local ( 2 H) of ~0.4 MHz was obtained using the same calculation, which agrees with the experimental T ( 2 H) of 0.27 MHz, because no p orbital is involved in the bonding of H/D ligands, thus the absence of local contribution to T ( 2 H).
Morton and Preston calculated isotropic hyperfine parameters for unit spin density (A) on an s orbital, and anisotropic hyperfine parameters (P) of unit spin density on a p, d, or f orbital for various elements: 14 where γ is the magnetogyric ratio of a given nucleus. By comparing our experimental Aiso and Tlocal values with those from Morton and Preston, it was possible to extract the electron density on each orbital: %s = Aiso/A × 100 and %p (or %d) = Tlocal/αP × 100, where α is the angular factor of the orbital (α = 2 /5 for p orbitals, and 2 /7 for d orbitals). The resulting electron densities are listed in Table 1 in the manuscript. The AFrame denotes the Euler angles, which define the orientation of the A tensor with respect to the g tensor in the convention used in EasySpin software. 9 Multiple quantum transitions were observed as a result of the large quadrupole moment of the 127 I nuclei and are plotted as red, yellow, green, and purple traces for ΔMI = ±2, ±3, ±4, and ±5 transitions, respectively. (Right) Echo-detected field-sweep D-band spectrum of 1-Y (black trace). The spectrum was fit to two paramagnetic components as indicated with red and yellow traces, the latter of which was assigned to an impurity constituting ~15% of the overall signal resulting possibly from air exposure of the sample. Arrows indicate the magnetic fields at which the EDNMR spectra (left) were measured.

Computational Methods and Analysis
For calculation of the IR spectra of 1-Y and 1-Y D , geometry optimization and calculation of the normal modes of vibration were performed in the gas-phase using unrestricted DFT within the Gaussian 09 rev. D package. 15 The single-crystal structure of 1-Y was used as a starting point, and all atomic positions were optimised simultaneously. The PBE density-functional was used in conjunction with Grimme's D3 dispersion correction, 16-20 the cc-pVTZ basis set was used for C, cc-pVDZ for H atoms, cc-pVTZ-PP for I atoms, 21,22 while the Stuttgart RSC 1997 effective core potential (ECP) was employed for the 28 core electrons of Y and the remaining valence electrons were described with the corresponding valence basis set. 23,24 Upon optimisation, the structure of 1-Y is modified only slightly, with an overall root-mean-squared difference of 0.198 Å across the entire structure (Table S8). The vibrational modes of the deuterated isotopologue 1-Y D were calculated by substituting the three hydride ions for deuterides (using the masses defined within the Gaussian code) and repeating the frequency analysis portion of the calculation in Gaussian using the same Hessian matrix as for 1-Y (Table S9). The theoretical IR spectra were generated by summing over the contributions of all modes { }, where each mode is represented by a Gaussian function centred at the DFT calculated mode frequency , with an area equal to the calculated linear absorption coefficient, and a full-width-at-half-maximum of 15 cm −1 .
The calculated IR spectra of 1-Y and 1-Y D are identical save for a handful of signals in the region 300-1500 cm −1 ( Figure S48) and agree well with the experimental spectra, which also only differ within this region ( Figure S11). The Vibrational Projection Analysis (ViPA) technique allows us to quantify the similarity between the normal modes of 1-Y and 1-Y D (Table S10). 25    For calculation of the UV-Vis spectra and hyperfine coupling of 1-Y (and 1-Y D ), and the exchange coupling of 1-Gd, the positions of the bridging hydride atoms were optimized from the respective X-Ray crystallographic structures using density-functional theory with ORCA v5.0.0. [26][27][28] In the case of 1-Gd, the Gd atoms were replaced with closed-shell Y to facilitate convergence of the SCF. The PBE0 density-functional was used in conjunction with the RI-COSX approximation, where 28-electron def2-ECP were employed for I and Y, and the def2-TZVP valence basis set was used for all atoms. 21,[29][30][31][32] For determination of the UV-Vis spectrum of 1-Y, the electronic structure was calculated using CASSCF methods in OpenMolcas v21.06, where the two-electron integrals were decomposed with a Cholesky decomposition (10 -8 threshold), and the basis sets were from the ANO-RCC library with VDZP quality for Y, I, H(hydride) and C(Cp), and with MB quality for all other atoms. [33][34][35] We employed the second order DKH relativistic decoupling. 36 State-averaged CASSCF (SA-CASSCF) calculations were performed for the 16 lowest doublet states in an active space of 1 electron in 16 orbitals, comprising orbitals with significant Y(4d) character. The ground state SOMO is a σ-like bonding orbital comprised of the dz 2 orbitals of the Y ions, while the space of Y(4d) excited states comprise various combinations with π-like and δ-like character (Table S11). Corrections for dynamic electron correlation were added using multiconfigurational pair-densityfunctional theory (MCPDFT) using the tPBE functional. 37 The UV-Vis spectrum was calculated on the basis of the SA-CASSCF transition intensities in the velocity gauge and the MCPDFT energies, and is in excellent agreement with the experiment ( Figure S49). The feature at ca. 13,200 cm −1 is the σ to σ* transition (state 1 to 4/5), and the feature at ca. 33,300 cm −1 is a σ to π transition (state 1 to 15).
For determination of the hyperfine coupling and spin density in 1-Y and 1-Y D , we employed ORCA v5.0.0 using the PBE density-functional with the ZORA relativistic Hamiltonian, the RI approximation, the ZORA-def2-TZVP basis set for all C and H atoms, and the old-ZORA-TZVP S62 basis set for Y and I. Isotropic and anisotropic hyperfine couplings were using the picture change corrections (Table S12) For determination of the exchange coupling in 1-Gd, we employed broken-symmetry DFT and SA-CASSCF-MCPDFT calculations. Broken-symmetry DFT calculations were performed with the Gaussian 09 rev. D package. 15 The B3LYP density-functional was used [38][39][40] in conjunction with the 6-31G* basis set for carbon and hydrogen atoms, 41 the 46 core electron Stuttgart-Dresden ECP and the corresponding double-zeta valence basis set for I, [42][43][44] while the 46 core electron Cundari Stevens Double Zeta ECP and the corresponding valence basis set was employed for Gd. 45 Calculations were performed for the high-spin, and four broken-symmetry solutions (Table S14 and Figure S51), leading to the exchange coupling values of J4f-σ = 176 cm −1 with JGd-Gd = −2.2 cm −1 obtained using the Noodleman method. 46 CASSCF-MCPDFT calculations were performed with OpenMolcas v21.06, 33 using a Cholesky decomposition for the two-electron integrals (10 −8 threshold), second order DKH relativistic decoupling, 36 and basis sets from the ANO-RCC library with VDZP quality for Gd, I, H(hydride) and C(Cp), and with MB quality for all other atoms. [34][35] The Gd atoms were first replaced with Lu (ANO-RCC-VDZP basis) and the σ-like SOMO was obtained using a CAS(1,1)SCF calculation ( Figure S50). The inactive orbitals were localized using the Pipek-Mezey method, 47 and the 4f orbitals for each Ln site were identified. CASSCF calculations were then performed for each Gd site separately, where the other two Gd sites remained defined as Lu, considering one S = 4 and one S = 3 root for an active space of 8 electrons in 8 orbitals, comprising the Gd 4f orbitals and the σ-like SOMO obtained from the Lu3 calculation which was not allowed to relax in the SCF procedure. Corrections for dynamic electron correlation were added using multiconfigurational pair-density-functional theory (MCPDFT) using the tPBE functional. 37 The energies show that J4f-σ is ferromagnetic and the average values for CASSCF and MCPDFT are J4f-σ = 333 cm −1 and 219 cm −1 , respectively (Table S15). Table S11. SA-CASSCF-MCPDFT energies (in cm −1 ) and natural orbitals for the low-lying excited states in 1-Y.  ↑ -↑ -↑ -↓ 3698 Figure S47. Illustration of the spin-spin coupling topology of the four spin centers in 1-Gd. The spin-spin exchange coupling values were obtained using the Noodleman method. 46 The calculated J4f-σ value agrees well with the value of J4f-σ = 168(1) cm −1 determined from experimental data as discussed in the main text. A small JGd-Gd = −2.2 cm −1 was also determined using DFT. Table S15. CASSCF and CASSCF-MCPDFT results for exchange coupling in 1-Gd. Figure S48. σ bonding SOMO of 1-Gd as determined from CAS(1,1)SCF calculations where Gd has been replaced with Lu.