Chemical tunnel-splitting-engineering in a dysprosium-based molecular nanomagnet

Total control over the electronic spin relaxation in molecular nanomagnets is the ultimate goal in the design of new molecules with evermore realizable applications in spin-based devices. For single-ion lanthanide systems, with strong spin–orbit coupling, the potential applications are linked to the energetic structure of the crystal field levels and quantum tunneling within the ground state. Structural engineering of the timescale of these tunneling events via appropriate design of crystal fields represents a fundamental challenge for the synthetic chemist, since tunnel splittings are expected to be suppressed by crystal field environments with sufficiently high-order symmetry. Here, we report the long missing study of the effect of a non-linear (C4) to pseudo-linear (D4d) change in crystal field symmetry in an otherwise chemically unaltered dysprosium complex. From a purely experimental study of crystal field levels and electronic spin dynamics at milliKelvin temperatures, we demonstrate the ensuing threefold reduction of the tunnel splitting.

χMT products and molar magnetization isotherms for 1Dy and 2Dy. The solid lines represent the global best fit as described in the main text. The experimental data were originally reported elsewhere 1 . The applied scaling factors (in order to account for error in the mass of the sample used) are 1.002 and 1.012 for 1Dy and 2Dy, respectively.   Figure 6 Selected part of the INS spectra for 2Dy (grey squares) and 2Y (orange squares) at T = 4 K recorded with an incident neutron wavelength of λi = 1.6 Å. The data were integrated over 1.2 Å -1 ≤ Q ≤ 2.5 Å -1 . The spectrum of 2Y can be well accounted for by a sum of two Gaussians (orange line) in the energy window of interest. As 2Dy and 2Y are isostructural, the spectral difference is assigned to the presence of a weak CF excitation in the spectrum of 2Dy. At λi = 1.6 Å, the Gaussian linewidths (full width at half maximum) of the lower-lying CF excitations (the ones from the ground doublet to the first and second excited doublets) were determined to be ~11.2 cm -1 (not shown). The spectrum of 2Dy in the energy transfers between 55 and 85 cm -1 was consequently described by a sum of three Gaussians; the two from the deconvolution of the 2Y spectrum, and a third with a linewidth of 11.2 cm -1 . In the fit, the linewidths of the Gaussians accounting for the two phonon modes were fixed to the values determined for 2Y, but their areas and positions were allowed to vary freely. The best fit is given by the red line, and the individual Gaussians are given as green dashed lines (phonons) or a solid blue line (CF excitation), respectively. The best fit affords a position of the CF excitation of 71(3) cm -1 . This value agrees reasonably well with the separation between the ground and third excited doublet of 67 cm -1 determined from the CF modelling. As discussed in the main text, this excitation is observable by INS due to the 25 % |±9/2 contribution to the wavefunction of the third excited doublet. The data were recorded on TOFTOF.

2Y 2Dy
Intensity / arb. units Figure 7 INS spectra for 2Dy (grey squares) and 2Y (orange squares) at T = 50 K and λi = 1.3 Å (a) and at T = 20 K and λi = 0.84 Å (b). The data were integrated over 3 Å -1 ≤ Q ≤ 5 Å -1 and 2.6 Å -1 ≤ Q ≤ 4.6 Å -1 , respectively. The blue lines give the INS spectra simulated using the CF parameters described in the main text and the relevant experimental parameters. The peak in the simulated spectrum corresponds to INS allowed CF excitation from the first excited doublet, |±13/2, to the highest lying one, |±15/2. The weak intensity associated with the "hot" nature (i.e. originating from an excited state) of the excitation combined with the loss of the relevant low Q-region, and the significant background associated with the incoherent scattering of the hydrogen nuclei contained in the sample, account for the lack of success in observing the excitation. The difficulties associated with observing CF excitations at very high energy transfers for hydrogen containing molecule-based materials have been reported elsewhere 2 . In both figures, the errors are less than the size of the symbols. The data were recorded on IN4C. Supplementary Figure 9 Schematics of the experimental setup for the cantilever torque magnetometry experiment along with the definition of the θ angle. Due to the macroscopic shape of the crystal (thin square-shaped plate) it was only possible to perform the rotation from the ab plane to the c axis. A simulation of the in-plane (the ab plane that is) rotation (Supplementary Fig. 10) demonstrates that its intensity is at least three orders of magnitude lower than the rotation experimentally performed. This justifies our choice to fix the crystallographic abc reference frame to be coincident with the magnetic xyz reference frame.
Supplementary Figure 10 Simulation of the in-plane rotation for 2Dy. Note that the signal is at least three orders of magnitude smaller than that observed for the rotation performed experimentally (i.e. the rotation from the ab plane to the c axis, see Supplementary Fig. 8).  Figure 11 Frequency dependence of the ac molar magnetic susceptibility for 1Dy in zero applied dc field measured at selected temperatures between 1.9 K and 5.5 K. Solid lines are best fits to a Cole-Cole function.  ac / Hz Figure 12 Frequency dependence of the ac molar magnetic susceptibility for 2Dy in zero applied dc field measured at selected temperatures between 1.9 K and 6.7 K. Solid lines are best fits to a Cole-Cole function.   Figure 14 Frequency dependence of the ac magnetic susceptibility for 1Dy in zero applied dc field measured at selected temperatures between 52 mK and 700 mK. Solid lines are best fits to a Cole-Cole function. Within the Debye model, the lattice contribution to the molar specific heat (C/R, R being the gas constant) is given by 3 C/R = (12π 4 r/5)(T/θD) 3 , where r is the number of atoms per formula unit, T is the temperature, and θD is the Debye temperature. As the expression only applies for T << θD, a fit of the molar specific heat was only carried out for T < 9 K. The so-obtained Debye temperatures are 212(1) K and 183(1) K for 1Y and 2Y, respectively. The best fits are given as solid lines. A plot of C/R vs. T 3 highlights the Debye contribution to the specific heat at low temperatures (c). Based on their isostructural nature 1 and the strong similarities in the phonon part of the INS spectra of the yttrium and dysprosium derivatives (see Supplementary Figs 3, 6-7), the Debye temperatures for 1Y and 2Y, determined above, should be reasonably applicable to the respective dysprosium derivatives. Within the Debye model, the mean phonon velocity, vm, in a given lattice can be calculated as 3 vm =(θDkB/h)(3NAρ/4πMw) -1/3 , where kB is the Boltzmann constant, h is the Planck constant, NA is Avogadro's constant, ρ is the crystal density, and Mw is the molar weight. From the Debye temperatures determined above, and the crystal densities determined from the X-ray single crystal structures reported elsewhere 1 , values of vm = 6.98(3)·10 3 m s -1 and vm = 6.55(3)·10 3 m s -1 are obtained for 1Dy and 2Dy, respectively. In order to correctly determine the dynamic crystal potential matrix element in the Orbach analysis of the spin-lattice relaxation rates of 1Y0.95Dy0.05 ( Supplementary Fig. 22), the mean phonon velocity of 1Y, also needs to be calculated. Using the appropriate molar weight and crystal density, a mean phonon velocity of vm = 6.98(3)·10 3 m s -1 is found for 1Y, the same as for 1Dy.  Supplementary Fig. 17. The field dependence of the spin-lattice relaxation rate can in the general case be described by the expression 4 T1 -1 (B,T) = B1/(1+B2B 2 ) + ATB n + k(T), where the first term accounts for tunneling, the second for direct relaxation between the two components of the ground doublet, and the last constant term absorbs any field independent relaxation (Orbach relaxation in this case). The best fit of the T1 -1 (B) data for 1Dy (a) to this expression with T = 5.0 K yields parameter values of B1 = 1.77(8)·10 3 s -1 , B2 = 4.6(8)·10 -6 mT -2 , A = 6.1(2)·10 -5 s -1 K -1 mT -2.6 , n = 2.6 (fixed), and k(T) = 3.22(4)·10 3 s -1 . For the standard direct process n = 4, but such a field dependence does not describe the high field data well. The reduced value of the exponent indicates that the direct relaxation process is phonon-bottlenecked, which is supported by the distinct increase in the value of α at high fields (b) 5 . The field dependence of the spin-lattice relaxation rate for 1Dy demonstrates that at zero field and temperatures as high as T = 5.0 K, tunneling within the ground state still influences the measured rates. Supplementary Figure 19 Frequency dependence of the ac molar magnetic susceptibility for 1Dy in an applied field of B = 120 mT measured at selected temperatures between 2.5 K and 6.1 K. Solid lines are best fits to a Cole-Cole function. The triangles give the spin-lattice relaxation rate for the concentrated specimen 1Dy, and the green line is the best fit of that data (equivalent to the best fit given in Fig. 4c of the main text).  Figure 24 Temperature dependence of the ac magnetic susceptibility in the mK regime, at the indicated frequencies, for 1Dy (a) and 2Dy (b). Temperature dependence of the isothermal (grey symbols) and adiabatic (open symbols) susceptibilities for 1Dy (c) and 2Dy (d) obtained from the Cole-Cole fitting of the frequency dependent data in Supplementary Figs 14 and 15. The appearance of maxima in the χ'(T) curves (a and b) signals the onset of long-range magnetic order. The setup used for the mK susceptibility measurements did not allow for a quantification of the susceptibility, thereby impeding an identification of the nature of the magnetic order through the comparison with a demagnetization factor. However, a Curie-Weiss analysis of the isothermal susceptibilities ( Supplementary Fig. 23) demonstrates the overall interactions between dysprosium ions to be antiferromagnetic in nature. The maxima in the χ''(T) curves suggest the presence of an uncompensated moment in the ordered phase, something that is incompatible with antiferromagnetic ordering given the symmetry of the crystals. However, as the spin dynamics are on the timescale of the experiment, these are in fact the origin of the maxima in the χ''(T) curves. As shown in panels c and d, the temperature dependence the isothermal susceptibility, obtained from Cole-Cole fitting of the frequency dependent data, mimics the temperature dependence of the χ'(T) curves while the adiabatic susceptibility is essentially temperature independent. As the dynamics are on the timescale of the experiment, the product 2πνacτ is not too far from unity, giving rise to a non-zero χ''(T) response following Supplementary Equation (3) (see Supplementary Note 1). As τ and α (1Dy: 0.25 ≤ α ≤ 0.39, 2Dy: 0.49 ≤ α ≤ 0.62) only show limited variation in the mK regime, the temperature dependence of the imaginary response must be dominated by that of the isothermal susceptibility. In this way, a maximum appears in the χ''(T) curves for 1Dy, while an onset of a peak is observed for 2Dy. Based on these considerations, we assign the appearance of maxima in the χ'(T) curves to the onset of long-range antiferromagnetic ordering. Note that all molecules are oriented with their fourfold axis (experimentally verified to be the easy axis of magnetization from torque magnetometry) along the crystallographic c axis. As the ground state g-tensors are highly anisotropic with g||>> g┴, the dipolar field arising from the neighboring dysprosium spins is expected to fulfill the condition Bdip,|| > Bdip,┴ at any given Dy site.

Exp. RT
Supplementary Table 1 |mJ contributions to the wavefunctions of the crystal field states within the ground Russell-Saunders multiplet for 1Dy. Only contributions of 1 % or larger have been included. The third column gives twice the expectation value of the projection of the total angular momentum operator along z (2Jz), calculated using the relation JzgJµB = (-∂E/∂B)µB, where gJ is the Landé g-value (gJ = 4/3 for Dy 3+ ) and -∂E/∂B is computed as part of the MagProp output. In the case of a diagonal CF Hamiltonian Jz = mJ. In order to make the dominantly contributing |mJ state easily identifiable, the value of 2Jz for each eigenstate is reported below. Supplementary Table 2 |mJ contributions to the wavefunctions of the crystal field states within the ground Russell-Saunders multiplet for 2Dy. Only contributions of 1 % or larger have been included. The third column gives twice the expectation value of the projection of the total angular momentum operator along z (2Jz), calculated using the relation JzgJµB = (-∂E/∂B)µB, where gJ is the Landé g-value (gJ = 4/3 for Dy 3+ ) and -∂E/∂B is computed as part of the MagProp output. In the case of a diagonal CF Hamiltonian Jz = mJ. In order to make the dominantly contributing |mJ state easily identifiable, the value of 2Jz for each eigenstate is reported below. Supplementary Note 1 All frequency dependent ac magnetic susceptibility data were fitted using a Cole-Cole function 6 in which the susceptibility is described by where νac is frequency of the oscillating magnetic field, χT the isothermal susceptibility, χS the adiabatic susceptibility, τ the relaxation time (i.e. the spin-lattice relaxation time T1 if the spin system is at thermal equilibrium with the phonon bath), and α the relaxation time distribution parameter. The latter parameter can take on values of 0 ≤ α < 1, where α = 0 corresponds to a single relaxation time. Supplementary Equation (1) can be rewritten into χ(νac) = χ'(νac) + iχ''(νac), where the real (in-phase) and imaginary (out-of-phase) components are given by where in our case, τ = T1.
Supplementary Note 2 In the expression used to describe the temperature dependence associated with the Orbach relaxation process (second term of Equation (2) in the main text), the parameter MO 2 appears. This parameter reflects the magnitude of the matrix elements of the dynamic crystal field potential relevant for the individual steps of the relaxation process. By combining the notation of Bierig, Weber, and Warshaw 7 and Larson and Jeffries 8 , MO 2 is defined as 2 2 ,  (4) should in principle contain an extra term with |d substituting for |c. However, as the terms are equal, the contribution can be accounted for by a factor of two 8 , which has been included in Equation (2) of the main text. In the phenomenological approach developed by Orbach, the thermal deformation of the crystal field is possibly of such low symmetry that all bk q might be non-zero 8 For 2Dy on the other hand, the ΔmJ = 12 difference between the |+13/2 state and the |-11/2 component of the ground doublet necessitates the inclusion of the 6 As outlined above, the g6 5 and g6 6 normalizing factors connecting the dynamic CF parameters to the diagonal static sixth order parameter is expected to be of the same order of magnitude but different. As the dynamic sixth order CF parameters furthermore are related to the determined values of B6 0 , the observed difference in the zero field values of MO 2 is the expected result.
Supplementary Note 3 Calculations of the tunnel splittings were carried out within the framework of Prokof'ev-Stamp theory 9-12 . In the low-temperature regime above TN, weak spin-spin interactions will give rise to a distribution of dynamic dipolar fields at each dysprosium site. The distribution width can be estimated from the previously determined Néel temperatures according to 11 For gJ = 4/3 and J = 15/2, values of dip = 4.8 mT and dip = 3.7 mT are calculated for 1Dy and 2Dy, respectively. The dipolar fields cause a splitting of the ground state of E = (dip 2 + T 2 ) 1/2 . dip is the dipolar energy bias, and T is the tunnel splitting in the ground state. The magnitude of the separation between the ground state and the first excited state, justifies an effective Jeff = ½ treatment in the temperature regime relevant for tunneling (at T = 5 K, the population of the first excited doublet is 0.1 ‰ and 0.2 ‰ for 1Dy and 2Dy, respectively). Our crystal field model affords the following components of the effective axial ground state g-tensors: 1Dy, g|| = 14.38, gc = 0.226; 2Dy, g|| = 14.53, g┴ = 0.081. By estimating the dipolar field experienced at a given site to be isotropic and estimated by the distribution width (i.e. Bdip,|| = Bdip,┴ ~ dip) 11 while the tunnel spitting is calculated as The values calculated for 1Dy are dip = 0.032 cm -1 and T = 5.1·10 -4 cm -1 , while dip = 0.025 cm -1 and T = 1.4·10 -4 cm -1 are obtained for 2Dy. The calculated ground state tunnel splitting for the lower symmetry derivative 1Dy is almost four times larger than the one in 2Dy.