Chemical control of spin–lattice relaxation to discover a room temperature molecular qubit

The second quantum revolution harnesses exquisite quantum control for a slate of diverse applications including sensing, communication, and computation. Of the many candidates for building quantum systems, molecules offer both tunability and specificity, but the principles to enable high temperature operation are not well established. Spin–lattice relaxation, represented by the time constant T1, is the primary factor dictating the high temperature performance of quantum bits (qubits), and serves as the upper limit on qubit coherence times (T2). For molecular qubits at elevated temperatures (>100 K), molecular vibrations facilitate rapid spin–lattice relaxation which limits T2 to well below operational minimums for certain quantum technologies. Here we identify the effects of controlling orbital angular momentum through metal coordination geometry and ligand rigidity via π-conjugation on T1 relaxation in three four-coordinate Cu2+S = ½ qubit candidates: bis(N,N′-dimethyl-4-amino-3-penten-2-imine) copper(ii) (Me2Nac)2 (1), bis(acetylacetone)ethylenediamine copper(ii) Cu(acacen) (2), and tetramethyltetraazaannulene copper(ii) Cu(tmtaa) (3). We obtain significant T1 improvement upon changing from tetrahedral to square planar geometries through changes in orbital angular momentum. T1 is further improved with greater π-conjugation in the ligand framework. Our electronic structure calculations reveal that the reduced motion of low energy vibrations in the primary coordination sphere slows relaxation and increases T1. These principles enable us to report a new molecular qubit candidate with room temperature T2 = 0.43 μs, and establishes guidelines for designing novel qubit candidates operating above 100 K.


References S85
S5 suppressing certain transitions and increasing the intensity of others. On those occasions, the glasses were remade prior to measurement. This preferred orientation was infrequently observed. For all three systems, T 1 was shorter in the glassy matrix than in the solid state at 5 K (5 ms in 1′ compared to 3 ms in 1′′, 14 ms in 2′ compared to 7 ms in 2′′, and 11 ms in 3′ compared to 1 ms in 3′′). However, these differences disappeared at higher temperatures. In the regime where relaxation was dominated by local mode relaxation, T 1 was largely independent of the matrix. This highlights an important facet of spin-lattice relaxation in S = 1 / 2 first row transition metal systems: low temperature dynamics are dominated by lattice phonons whereas high temperature dynamics are dominated by intrinsic molecular properties. At low temperatures, therefore, the unique phonon densities of different matrices can have large impacts on relaxation behavior. At high temperatures, local mode mediated relaxation dominates, and dynamics can be largely considered matrix agnostic. T m increases in OTP solutions due to the low nuclear spin density of OTP. 4 This effect is most strongly observed in 1, where removing the methyl group rich environment of 4 enhances T m by nearly 2 orders of magnitude (at 5 K, T m = 0.32 μs in 1′ compared to T m = 14 μs in 1′′).
We observe the same trends observed in the relaxation dynamics of 1′-3′ in 1′′-3′′ (Table S7). However, there are some slight changes in magnitude of the relaxation parameters. A Dir increases in OTP, which is expected from the shorter T 1 in the glass. While we observe a change in B Ram , it is important to note that B Ram can be selected somewhat arbitrarily for a given value of Θ D . 5 Since the Θ D of the OTP glass was assigned an arbitrary value of 65 K for all three fits, whereas Θ D in the solid state systems was a fit parameter unique to each system, changes in B Ram could be attributed to changes in Θ D that we cannot model. Deviations in C Loc and Δ Loc are small and frequently within the error of the fit, further reinforcing that local mode relaxation is matrix agnostic.
Electron Paramagnetic Resonance Measurements. Samples were prepared for analysis by a solid-state dilution in their respective diamagnetic analogues. Solid-state dilutions were prepared in a 1:100 (1%) ratio to suppress the influence of intermolecular electronic spin interactions on T m and T 1 . All samples were loaded into 4 mm OD quartz tubes (Wilmad 707-SQ-250M), restrained with eicosane, and flame-sealed under high vacuum. Solution phase measurements were prepared in a 1% by weight dilution with OTP. Though this leads to varying spin concentrations of spin/unit volume, the variation in terms is small (between 1 and 4% spin/volume). The effects of these concentration differences are minimal, and our T m times are more heavily influenced by the nuclear spin concentration of the matrix (as discussed above). Additionally, the effects of cross relaxation are naturally suppressed at the higher temperatures of interest in this study (70 K). 6 Spectroscopic data were obtained on either a Bruker E580 X/Q-band spectrometer or a Bruker E560 X-band spectrometer, both with a split ring resonator (ER4118X-MS5) and a 1 kW TWT amplifier (Applied Systems Engineering). Prior to all pulsed measurements, the resonator was over-coupled to minimize ringdown following application of the microwave pulses. Temperature was controlled with an Oxford Instruments CF935 helium cryostat and an Oxford Instruments ITC503 temperature controller. All EPR data were processed by a combination of XEpr (data collection), 7 Matlab (phasing and normalization of EPR data), 8 Easyspin (modeling and simulation of cw-EPR spectra), 9 and Origin (fitting of decay curves). 10 Absolute intensities of the cw-EPR spectra were normalized between 0 and 1 then, simulated using the pepper function in EasySpin.
Spin-lattice relaxation times were obtained using saturation recovery sequences at the highest intensity peaks in the echo-detected spectrum. These sequences achieved saturation by applying a picket-fence S6 saturation sequence of twenty consecutive 12 ns pulses. Following a delay time T beginning at 100 ns, Hahn echo detection was used to monitor the recovery from saturation with π / 2 and π pulses of 16 and 32 ns, respectively. Four-step phase cycling was used on these measurements. All data was phased by maximization of the sum of the square of the data in the real component of the spectrum. For 2′-3′ and 1′′− 3′′, data was collected with logarithmic spacing, allowing for data to quickly be taken quickly across orders of magnitude change in time. The T 1 of 1′ was collected using a slightly different detection scheme wherein data was collected with linear spacing, allowing us to capture more points in the regions of largest change. The best fit was obtained using a monoexponential fitting function 11 using the following equation, . The fit to low temperature data could be improved by including a stretch factor (σ) to account for the influence of multiple processes with slightly different T 1 constants using the equation . This distribution of T 1 relaxation times has been attributed to a handful of factors, most likely in this report are inequivalences between magnetic cites, causing slight deviation in relaxation rates. 12 A Stretch factor was often required to adequately fit low temperature data. At higher temperatures where single phonon relaxation becomes less prominent, these inequivalences matter less and the stretch factor trends towards 1 as shown in Tables S14-S19.
The temperature dependence of T 1 was fit to account for the direct process, Raman process, and local modes using the following expression, Where T is temperature, A Dir is the coefficient for the direct process, B Ram is the coefficient for the Raman process, J 8 is the transport integral, 13 Θ D is the Debye temperature, C Loc is the coefficient for local modes, and  loc is the energy of the active local modes of vibration. See below for further discussion of the fitting process of the temperature dependence of T 1 across 1′-3′ and 1′′−3′′. Phase memory times (T m ) were obtained using a two-pulse Hahn echo sequence, π / 2 -τ-π-τ-echo, where τ is the time delay between pulses, and π / 2 or π denote microwave pulses, 16 ns and 32 ns, respectively. Starting delay times were selected to minimize the effects of any observed ringdown. These delays ranged from 100 to 200 ns for all complexes. Four step phase cycling was employed for these experiments as well.
Most systems were modeled using a simple monoexponential decay function . Some systems also exhibited electron spin echo envelope modulation (ESEEM) of the echo intensity, primarily from weakly coupled 1 H nuclei. Due to the inherent errors in adding the necessary fit parameters to the relaxation equation to model ESEEM, we chose to ignore the ESEEM in our fits. Nutation experiments were performed by applying a tipping pulse (τ p ) to generate a superposition state, followed by followed by a Hahn-echo detection sequence (τ p -T-π / 2 -τ-π-τ-echo). For constant τ and microwave power, the echo intensity varies with the length of the tipping pulse. As τ p is incremented, the spin samples all superposition states, resulting in a sinusoidal oscillation of echo intensity. Plotting the oscillation in echo intensity versus the tipping pulse length gives a decaying sinusoidal curve, known as a Rabi oscillation -the hallmark of qubit viability. 14,15 For our experiments, τ p was incremented in 2 ns steps S7 with T = 340 ns, and π / 2 and π pulse lengths of 16 ns and 32 ns respectively. For the room temperature measurements of 3′, a large baseline needed to be subtracted from the data. This baseline was only present for these high temperature measurements, suggesting it results from the weak echo intensity.
The Fourier transform of nutation data yields plots of the frequency of oscillation in the experiment. At low temperatures, we observe a secondary peak at approximately 14 MHz, the Larmor frequency of the nearby proton nuclei, as a result of the Hartman-Hahn effect. 16 To demonstrate room temperature qubit manipulation of 3′, we measured the rate of Rabi oscillations (known as the Rabi rate) at three microwave powers. The relative microwave power was calculated as: where A is the microwave power attenuation and Z is the lowest microwave power attenuation in the set of three molecules examined in the experiment (7.2 dB). Although the trend between Rabi rate and relative microwave power attenuation is expected to be linear, we observe a small deviation from linearity (R 2 = 0.988). We attribute this to the weak echo intensity of 3′ at low microwave powers. Because the echo intensity is so weak, we are only able to observe one period of the oscillation, leading to error in our estimation of its Rabi rate.
X-ray Structure Determination. All diffraction data were collected in the X-ray crystallography lab of the Integrated Molecular Structure Education and Research Center (IMSERC) at Northwestern University. Single crystals of 1, 2, 4, and 5 suitable for X-ray diffraction analysis were coated in Paratone N oil and mounted on a MiTeGen MicroLoopTM. Crystallographic data for 1 and 2 were collected on a Rigaku diffractometer equipped with a MoKα sealed tube diffraction source with a graphite monochromator, and a Bruker APEX II detector. Data for 4 and 5 were acquired on a Rigaku XtaLab Synergy, single source at home/near, HyPix diffractometer equipped with a MoKα PhotonJet X-ray source. All datasets were collected at 100 K. Raw data were integrated and corrected for Lorentz and polarization effects with SAINT v8.27B.10 (for 1 and 2) or CrysAlisPro 1.171.40.53 (for 4 and 5) 17 Absorption corrections were applied using SADABS. 18 Space group assignments were determined by examination of systematic absences, Estatistics, and successive refinement of the structures. Structures were solved using direct methods in SHELXT and further refined with SHELXL-2013 19 operated with the OLEX2 interface. 20 Hydrogen atoms were placed in ideal positions and refined using a riding model for all structures. Structural parameters for 1-2 were obtained from single crystal X-ray diffraction experiments. (Table S1-2).Structural parameters for 3 were obtained from the literature reported crystal structure. 21 To ensure crystallographic similarity, powder X-ray diffraction samples of 3 were compared to simulated powder X-ray diffraction patterns obtained from Mercury with a full-width half-max of 0.5° 2θ. A well-ground powder sample of 3 was sandwiched between pieces of Kapton tape. Powder X-ray diffraction patterns were collected on a STOE STADI MP diffractometer equipped with CuKα radiation. The provided crystal structure of 1 shows a positional disorder of the complex in the unit cell with approximately 90% of the complexes of one orientation, and 10% with the antiparallel orientation. Solving the structure without the inclusion of this second orientation leads to a significant amount of unmodeled electron density, but because of the low occupation, there is significant disorder in modeling this orientation, resulting in an A-level and B-level alert. In 2, there is a difficult to model, disordered, partially occupied water molecule in the crystal structure. The exact orientation of the oxygen atom causes a B-level alert, likely due to occupational disorder. In 2, additional alerts were caused by the size of the crystal. Due to the thinness of the crystal, we decided not to cut it in order to maximize our sample exposure to X-rays. This did however lead to a larger than ideal crystal, causing an A-level alert.
The reported structure of 3 was taken at room temperature (298 K) and is being compared to structures at 100 K. Thought this temperature difference does make exact comparison of bond lengths difficult -our assignment of relative τ 4 ' values are validated by comparison of the computed relaxed geometries (Table  S20).
To confirm our determined crystal structure, we performed powder X-ray diffraction experiments on 1-3 at room temperature with CuKα (1.54065 Å) radiation (Figure S5-S7). The resulting diffraction patterns were compared to the predicted spectra from their single crystal structures. There is a slight shift in peak 2θ present in 1-2, resulting from comparing powder patterns obtained at 300 K to crystal structures obtained at 100 K.
Other Physical Measurements. Infrared spectra were recorded on a Bruker Alpha FTIR spectrometer equipped with an attenuated total reflectance accessory and diamond anvil. Electrospray ionization mass spectrometry measurements were performed on MeOH/DCM solutions of 1-6 with Bruker AmaZon SL ESI-Ion Trap Mass Spectrometers at the IMSERC facility of Northwestern. UV-Vis spectra were collected on toluene solutions of 1-6 with a Varian Cary 5000 spectrophotometer.
Discussion of sources of error in fitting variable temperature T 1 . In modeling these data, we focused on the direct process, the Raman process, and the local mode process. As all complexes are monometallic and S = ½, an Orbach mechanism, which necessitates accessible low-lying spin states only available in higher spin systems or exchange coupled systems with multiple spin centers, is physically unreasonable and would lead to overparameterization. 22,23 Because the model used to fit the data (The Debye model) was derived for extended solid systems of much higher crystallographic symmetry (the highest symmetry examined in this report is orthorhombic whereas the Debye model was developed for cubic systems) there is a large degree of inherent error associated with the fit parameters for these systems. 23,24 The Debye temperature obtained from these fits is not an analogous parameter to the Debye temperature of ideal solids. In these systems, it is far better understood as metric representing the phonon energy of an individual lattice. Additionally, the Debye temperature and Raman coefficient B Ram are heavily interdependent. It has been demonstrated that a change in B Ram of approximately one order of magnitude can be offset by a change in the Debye temperature by approximately 30%. 5 This is typically corrected for by limiting the range of acceptable Debye temperatures to a range typically associated with molecular solids (between 50 K -and 120 K). 22,24 We use this to our advantage in fitting the OTP data to hold lattice phonon contributions constant. We first fit the data for 3′′ with no parameters held constant, and with the aforementioned constraint on Θ D . This led to a fit Θ D of 65 K. We then used this as a representative Θ D for the OTP matrix, refit the data for 3′′ with Θ D now held constant at 65 K, as well as fitting the T 1 data for 1′′-2′′ as well.
Additionally, significant errors are associated with the direct process, because this mechanism is operative at the lowest temperatures of measurements and is typically fit to only the few data points below 10 K. The Δ Loc parameter heavily discussed in the main text represents some weighted average of all local modes contributing to spin-lattice relaxation. By the strictest definition, the term represents a single mode in the system. However, as we and others demonstrate, modeling these systems with only a single impactful local S9 mode is not accurate. 25 Technically, one could add additional terms -one for each mode being modeled. 26 However, adding additional terms to the Debye model leads to problematic overparameterization. 27 Therefore, convention is to fit relaxation to a single Δ loc term, with the caveat that it is likely not representative of a true mode in the system, and instead gives information about the average of all modes in the system. As we demonstrate in the main text, this method does give partial information about vibrational modes responsible for relaxation, but does not paint a complete picture of the local modes process in molecular qubit candidates.
Modeling of cw-EPR Spectra cw-EPR spectra were modeled through the program EasySpin, run through MatLab. 8 In all systems, anisotropy in the g-tensors and hyperfine interaction (A-tensors) were considered. The hyperfine interaction is strongly distance dependence, and is generally only significant for the primary coordination sphere in transition metal systems. Therefore, the only atoms treated as having measurable hyperfine interactions in our model were the metal center and the nuclear spin bearing atoms of the primary coordination sphere (nitrogen in all systems). In an attempt to further reduce spin-spin broadening, we performed cw-EPR on 0.1 % by weight solutions of 1-3 in OTP (1′′′-3′′′). Linewidths remained largely consistent with the higher electronic spin concentration 1′-3′, suggesting cw-EPR linewidths are not limited by interelectronic interactions, and are instead limited by interactions with nearby nuclear spins. Both 1 and 3 were modeled with an anisotropic line broadening term. Easyspin supports modeling anisotropic strain with several different methods. To reduce the parameter space, we chose to model the systems with only a single source of line broadening. It should be noted, however, that likely a combination of many sources of anisotropy contribute to line broadening in each system. For 1′ and 1′′′, the most success was found using an "H-strain" term, an anisotropic smearing of magnetic parameters within the system causing broadening of the observed transitions. The choice to model with H-strain was made as there is no clear chemical origin of the anisotropy.
3′ was best modeled with a "g-strain" parameter, which models an anisotropic distribution of g tensors existing in the crystal causing an anisotropic broadening of transitions. In the case of 3′, the use of g-strain makes particular sense as, in the crystal structure of both 3 and 6, the molecule occupies two unique positions in the solid-state crystal. In 3′, 3 can occupy either of the two sites, and therefore has a cw-EPR spectra composited of 2 different molecules, with similar but distinct magnetic parameters. It is very likely that an even more accurate simulation could be made from specifically modeling two unique, uncoupled copper sites in the system. The existence of this broadening is highlighted by the significant decrease in anisotropic broadening in 3′′′, wherein the large "smear" of transitions around 340 mT is resolved into several distinguishable peaks. In 3′′′, the system could now be modeled with an isotropic Lorentzian.
Although the crystal structure of 2 also has two unique sites in the crystal, 2′ is well modeled with an isotropic Lorentzian. This is best understood as an effect of the different matrices the spins are in. Both 3 and 6 crystalize into the P2 1 /c space group with two unique molecule sites per unit cell. 21,28 Although 2 and 5 both crystalize into the C2/c space group, the unit cell of 5 has one molecular site in the cell.. When 2 is diluted into 5 to make 2′, it substitutes into the crystal structure of 5, and there is only 1 magnetic site the molecule can occupy. Therefore, strain effects are not observed, and we obtained the well resolved cw-EPR spectra reported in the main text (Figure 2). 2′′′ was similarly modeled with a simple isotropic Lorentzian.

Notes on Comparing Geometries of Different Complexes
The dihedral angle cited in the main text and in figures was defined as follows: A single subunit is defined as a 5-carbon chain with a main group donor atom located on the 2 and 4 carbons. In all three complexes, a plane was defined as the plane created by the two donor atoms on a NacNac or NacAc subunit and the metal center. Since each complex contained two subunits, the dihedral angle discussed in the paper is the angle between these two intersecting planes. As one would expect, it is large for a tetrahedral or pseudo tetrahedral molecule, and small for a square planar molecule.
There exist 2 common parameters which compare the geometry of four coordinate molecules to rank how similar they are to those geometric extremes -τ 4 and τ 4 ', defined as follows: Where α and β are the two largest valence angles in the complex. 29,30 Although the parameters are similar at the extreme (0 being perfectly square planar, 1 being perfectly tetrahedral), τ 4 ' uniquely identifies the α and β angles of the complex, and therefore better highlights the unique distortions between structures.
Computational Details Vibrational modes, stiffness, and spin density for each of the molecules was calculated using the Γ-point version of VASP [31][32][33] with the Perdew-Burke-Ernzerhof (PBE) functional. A plane wave cutoff of 600 eV, an augmentation grid plane wave cutoff of 2400 eV, and an energy convergence of 10 -8 eV was used for all VASP calculations. Molecular structures were extracted from experimental crystal structures, giving a gas phase approximation. The structures were relaxed with a force convergence of 5x10 -5 meV/Å and a minimum of 10 Å gap between the molecule and its periodic images was used. The vibrational modes were calculated from a dynamical matrix populated using finite displacements of 0.01 Å of each atom. The vibrational density of states was calculated via a gaussian folding method with a width of 4 cm -1 . 34 The intensities of each mode were normalized to one. Eigenvalues of each mode were projected onto individual atoms and then summed by atom type to give the full atom projected density of states.
Local stiffness of the Cu atoms was found by fitting the local potential energy , with the equation: (1) Where is reference potential energy, is the finite displacement vector from the equilibrium position, 0 ∆ and is the matrix of second order polynomial terms approximating the Hessian matrix. Finite displacements were evaluated on a 5 x 5 x 5 grid around the equilibrium position in increments of 0.01 Å. Spin density is shown with an isosurface containing 70 percent of the density.
The structures were re-relaxed with ORCA 35 version 4.0.1.2 using the B3LYP functional [36][37][38][39] and the def2-TZVP 40 basis set with the def2-/JK 41 auxiliary basis set. Hyperfine coupling and g-tensor calculations were performed on these re-relaxed structures using the ORCA program. All g-tensor calculations used the same functional and basis set as the relaxations. For hyperfine coupling the functional and basis set remained the same as other ORCA calculations with the exception that the triply polarized basis set CP(PPP) 12 was used for the Cu atom. For all ORCA calculations, tight SCF convergence criteria (giving an energy change of 2.7211x10 -7 eV) and DFT integration grid 6 (Lebedev 590 points for the angular grid and S11 the IntAcc parameter set to 5.34 for the radial grid) were used. A tighter DFT integration, grid 7, was used around the Cu atom (Lebedev 770 points for the angular grid and the IntAcc parameter set to 5.67 for the radial grid). These setting allowed for calculation of g-tensor values that were converged to 1x10 -4 . The g-tensor and hyperfine couplings to the vibrational modes were calculated for each normal mode with five amplitudes, four between 0.02 and 0.08 Å. One additional amplitude at -0.08 Å was also computed to indicate if the changes in the g-tensor with respect to mode amplitude were even or odd for a given mode. All g-tensor data for distorted structures was fit to a linear model and a quadratic model to determine the best fit for each mode. For the linear model the data was fit to Equation 2, where g eq is the equilibrium gtensor value, d amplitude is the amplitude of the distortion applied to the equilibrium structure, and g dis is the gtensor value for the distorted structure. (3) The calculated g-tensors for each distorted structure showed that modes with significant coupling, greater than 10 -3 Å -1 , were well described by a linear approximation with most modes having an r 2 value above 0.95 (Table 27-S29). Several of these strongly coupled modes for 1 show a "V-shaped" behavior and are identified in Table S27. Modes that were well described by a quadratic approximation all have coupling constants below 10 -3 Å -1 . Occasionally a mode could not be fit to a quadratic approximation due to the value being constant with respect to displacement. Displacements where this occurs are notated with "nan" in Table S27-S44. The linear approximation of the first derivative of the g-tensor with respect to the phonon modes, was then used the calculate a weighted spin phonon coupling (Equation 4) Where ω, is the frequency of the mode and k b is Boltzmann's constant. Temperature, T, ranges from 0-300 K and Q ranges from 0 to n where n is the number of real modes in the frequency range 0-500 cm -1 . This range was chosen as representative of all modes with a realistic thermal occupancy.
To analyze the perturbation of the primary coordination sphere during relaxation, we applied a bond distance that computes the average distance between the copper cation (C) and every nitrogen donor (N) with respect to each axis (Equation 5).
This metric was calculated for every distorted structure and then fit to a linear model with the distortion amplitude to determine the change in the bond length metric with respect to distortion amplitude. Larger d Cu-N values indicate more distortion in the primary coordination sphere along a vibrational mode. For 2, this metric was altered slightly to include the perturbation in the donor oxygen atoms as well. The description and calculation are identical, treating oxygen displacement identically to nitrogen displacement.
To identify the unique environment, the parameters was termed d Cu-E .
A simple thermal model based off the work of Kazmierczak et al 42 was used to calculate T 1 for all three molecules. The model weights the SPC coefficients by their Bose-Einstein occupancy and then is normalized with experimental T 1 data for Cu(tmtaa) in o-Terphenyl (eq. 6). A constant normalization factor (a) was determined from the data at 160K.

Comparison of PDOS with Experimental Raman and IR Spectra:
The gas phase approximation was used for phonon calculations. This has been demonstrated to cause discrepancies at very low frequencies due to the lack of lattice modes present in the approximation. 27 Direct comparison of predicted vibrational behavior to experimental results is difficult. While calculations necessitate ignoring the local matrix environment to prevent the exponential scaling of computational difficulty, these factors necessarily impact vibrational modes, both in energy and in oscillator strength. As such, we wish to stress that our results from DFT calculations cannot be 1:1 mapped onto experimental results in real matrices. For example, the vibrational modes predicted at 129 cm -1 for 2 that is heavily discussed in the main text, may be shifted by tens of wavenumbers in energy in a real matrix, and these shifts could also shift the predicted spin-phonon coupling of those vibrational terms. As such, we only use our model to guide our understanding of experimental results, and stress that a direct comparison of experimental and theoretical vibrational energies is non-trivial.
Despite these approximations, most Raman modes match closely with a computed mode, giving an overall low root mean square error of 19.73, 17.13, and 11.74 cm -1 for 1-3 respectively. This error is similar to that found in other four coordinate copper complexes. 43 S13 where n is the number of reflections and p is the total number of parameters refined where n is the number of reflections and p is the total number of parameters refined where n is the number of reflections and p is the total number of parameters refined where n is the number of reflections and p is the total number of parameters refined   Figure S1 | Compliance of the copper metal center for 1-3. A larger value indicates that there is less resistance to motion of the copper along that axis. Though this metric does not fully capture the rigidity of ligand π-bonds, it does demonstrate how tightly bound the spin bearing metal center is in the binding pocket.   Figure S2 | Energy of the spin-up and spin-down molecular orbitals derived from the copper d-orbitals in 1 (in eV), and the percent character of the molecular orbital comprised of the atomic d-orbital. Numbers index molecular orbitals in ascending energy. The HOMO orbital is assigned to 0 eV. The displayed orbitals are the five orbitals with the highest %d-orbital character.
S21    Figure S4 | Energy of the spin-up and spin-down molecular orbitals derived from the copper d-orbitals in 3 (in eV) and the percent character of the molecular orbital comprised of the atomic d-orbital. Numbers index molecular orbitals in ascending energy. The HOMO orbital is assigned to 0 eV. The displayed orbitals are the five orbitals with the highest %d-orbital character.    N). As the antibonding character of the orbital the electron is in increases, the spin becomes increasing localized on the ligand, represented by a decrease in ρ Cu and an increase in ρ E,Avg . The decreased spin density on the ligand in 2 is a result of weaker resonance in the ligand. In viewing the ligands of 1−3, the acacen ligand in 2 has weaker resonance contributions from structures which place increased spin density on carbon, causing a comparatively greater ρ E,Avg than would otherwise be expected.                 Figure S5 | PXRD comparison between ground crystals of 1 and the calculated pattern from our obtained crystal structure. The experimental pattern was collected at 298 K, and the simulation is calculated from the single crystal data taken at 100 K. Powder diffraction data were collected with pure CuKα (1.54065 Å) radiation.
S52 Figure S6 | PXRD comparison between ground crystals of 2 and the calculated pattern from our obtained crystal structure. The experimental pattern was collected at 298 K, and the simulation is calculated from the single crystal data taken at 100 K. Powder diffraction data were collected with pure CuKα (1.54065 Å) radiation.
S53 Figure S7 | PXRD comparison between ground up crystals of 3 with the literature reported structure. 21 Both data sets are at 298 K. Powder diffraction were was collected with pure CuKα (1.54065 Å) radiation.

Figure S8
| A simplified view of the local environment of Cu(Me 2 Nac) 2 (1′) in the solid state. This model shows all hydrogen atoms (white) within a 10 Å distance (represented by the purple sphere) from a copper metal center (orange). We attribute the unexpectedly short T m for 1′ to this large proton density of more than 100 hydrogen atoms, many on free-rotating methyl groups.         S78 Figure S32 | Comparison of the experimental Raman spectroscopy normalized intensities, experimental FTIR spectroscopy percent transmission, and computed normalized phonon density of states for 1. FTIR and Raman spectra were taken on powder samples at ambient temperature.
S79 Figure S33 | Comparison of the experimental Raman spectroscopy normalized intensities, experimental FTIR spectroscopy percent transmission, and computed normalized phonon density of states for 2. FTIR and Raman spectra were taken on powder samples at ambient temperature.
S80 Figure S34 | Comparison of the experimental Raman spectroscopy normalized intensities, experimental FTIR spectroscopy percent transmission, and computed normalized phonon density of states for 3. FTIR and Raman spectra were taken on powder samples at ambient temperature.  16 Additional low frequency, weak intensity peaks are artifacts of data processing. Figure S38 | Plots of spin densities from VASP calculations at the 70% occupancy level for 1−3, with the spin densities located on the copper atom (ρ Cu ) and the average located on each donor atom (ρ E,Avg ). There is an increased spin density on donor atoms and on the copper center in 2 relatives to 3 due to the strong resonance effects present in 3 allowing for more delocalization of the spin onto the carbon backbone of the ligand. This places more spin density on the direct donor atoms in 2. A weaker, though still present, resonance effect in 1 is also at play, explaining the comparatively similar ρ Cu between 1 and 2. Strong interactions with the ligand, as well this resonant effect, explain the comparatively delocalized spins in this report relative to others in the literature. 23