Magnetization relaxation in the single-ion magnet DySc2N@C80: quantum tunneling, magnetic dilution, and unconventional temperature dependence

Quantum tunneling and relaxation of magnetization in single molecule magnet DySc2N@C80 is thoroughly studied as a function of magnetic dilution, temperature, and magnetic field.


Adsorption of DySc2N@C80 by the MOF DUT-51(Zr)
DySc2N@C80 in toluene was prepared for the DUT-51(Zr) adsorption experiment. The amount of the DySc2N@C80 was estimated with UV-vis-NIR absorption spectroscopy according to the Lambert-Beer's Law: A=lg(I0/I)=L•c•ε, where A is absorbance, L is the optical length (1 cm), C is concentration, and ε is molar extinction coefficient at 420 nm, 15.8 x 10 3 M -1 cm -1 , taken from the value of Sc3N@C80 based on their similar structure. 9 The concentration of DySc2N@C80 during the adsorption was monitored by UV-vis-NIR spectroscopy. The results are shown in Figure S5. After a fast drop of the fullerene concentration in solution caused by adsorption into the voids of DUT-51(Zr) during the first 24 hours, the steady concentration was achieved after ca 100 hours, as shown in Figure S6. The fullerene solution over the DUT-51(Zr) was refreshed four times reaching in the end the fullerene/DUT-51(Zr) mass ratio of ca. 1:100 (0.18 mg of DySc2N@C80 was adsorbed by 19.6 mg DUT-51(Zr)).   Figure S6. The concentration evolution of DySc2N@C80 for the adsorption presented with the absorption of the solution at 420 nm. The data was extracted from Figure S5. Absorbance at 420 nm (a.u.)

Time (min)
Absorbance at 420 nm S10 Dilution in polystyrene Figure S7. Magnetization curves (T = 1.8 K) and microscope images of DySc2N@C80 diluted with polystyrene and drop casted from different solvents: toluene (black), CS2 (red). Crystals of fullerene can be clearly seen in the film drop-casted from toluene, whereas the film obtained from CS2 is much more homogeneous.
Homogeneity of the has a strong influence on magnetic hysteresis. . Normalized magnetization curves of DySc2N@C80 powder (black) and single-crystal of DySc2N@C80•Ni II (OEP). Initially the crystal was aligned with the long side of the block along the direction of the magnetic magnetic field. Light green curve is the magnetic curve for this orientation of the crystal (cooled-down in zero field). When the crystal is cooled down in the presence of magnetic field, it rotated by ca 45⁰. Magnetization curve of the magnetically aligned crystal is shown in olive. Also shown for comparison is the magnetization curve for the non-aligned crystal, which magnetic field axis is scaled by a factor of 0.73, which is equivalent to the rotation of the magnetic field direction by 43⁰. The perfect match between magnetization curves of the magnetically-aligned crystal and of the scaled curve of the nonaligned crystal can be seen.

Determination or magnetic relaxation times with DC magnetometry
The two most common macroscopic approaches for determination of the system's characteristic times are: DC relaxation measurements (above 100 s, relaxation curves are recorded) and AC relaxation measurements (below 10 s, searching for the resonance behavior of the sample). The DC measurement of the magnetization decay curve is relatively simple and doesn't require large sample amounts (5-10 % compared to AC). However, difficulties arise during fitting of the experimental data as single-molecule magnets tend to exhibit a time-dependent decay rate. One of the reasons is arising from the evolution of internal dipolar fields in the sample during relaxation. Consequently a single exponential function often fails to describe the system's behavior. In rare cases double exponential description is sufficient 10 . However, in a general case the decay curve consists of an infinite number of exponentials. A characteristic value for the relaxation time distributions has to be derived. It becomes possible with a stretched exponential: Where and 0 are the equilibrium and initial magnetizations, respectively, 1 is a characteristic relaxation time and is an additional parameter that corresponds to the time-dependent decay rate −1 ~ −1 with = (0; 1). In the extreme case of = 1 one obtains a single exponential.
The stretched exponential can be represented in the form that corresponds to the continuous sum of exponential decays: Where ( , ) is a normalized probability distribution 11 .
For the first time the stretched exponential was proposed by Rudolf Kohlrausch in 1854 to describe the relaxation of charge from a glass Leyden jar (the precursor of modern capacitors). This function is commonly used in the studies of glassy materials, were heterogenic processes are expected. For example, the stretched exponential relaxation can be found in simulations of the 3D spin glass with large enough number of spins 12 . In the field of single-molecule magnets a time-dependent decay rate was found already for Mn12ac, the first molecular magnet, consequently a stretched exponential was used 13 . The detailed comparison of the most common characteristic time derivation strategies is presented below. The two extreme cases of relaxation experiments should be taken into account. The first one corresponds to short relaxation times, which are close to the limit of the DC measurement technique. The limitation arises from the characteristics for the superconducting magnet, used in the experiment. Prior to recording the relaxation curve one has to magnetize the sample in the saturation field and then sweep the field down and stabilize in at the desired value. The common change of the magnetic is about 2-3 Tesla. For the modern QD MPMS3 magnetometer the delay between the start of sweeping and the first measured point is between 15 and 30 seconds. Thus, samples with fast relaxation undergo a considerable demagnetization already before the actual measurement starts.
The consequences are different for a constant and a time-dependent decay rates. While for the first case a single exponential fitting of the whole decay curve or a part of it would provide the same result, for the second case a curve that misses the initial part would provide overestimated values due to a biased representation of the real time distribution. Considerable deviations are found for the times below 100 s (QD MPMS3), therefore such results are excluded.
On Figure S10 an example of a relatively fast decay is shown. For such cases the measurement over more than 10 characteristic times is affordable, which provides an easier estimation of the equilibrium magnetization. This is crucial for the correct fit as the result exhibits a strong dependence on the value of Meq. As a consequence the fitting algorithm has to be chosen carefully. By default most of the routines are dealing with minimization of the sum of squared residuals (the differences between a data points and the S15 corresponding model estimates). This leads to unequal treatment of the "first" and the "last" points of the curves, as their values can differ by several order of magnitude, still having the same instrumental error.
To ensure that the average residual can't surpass the small values of the "last" points, the "normalized" residuals have to be introduced, calculated as the ratio between the residual and the corresponding value. Effectively this forces the fitting curve to follow the tail of the experimental curve. According to Figure S10 the stretched exponential describes the magnetization decay in DySc2N@C80 single crystal significantly better than a single exponential, and the "normalized" residuals approach increases the fit quality for the last points.
An alternative way of treating the time-dependent relaxation rate can be found in literature. The evolution of internal dipolar fields in the sample during relaxation can lead to a square-root time dependence of the signal 14 . However, this special type of relaxation is expected only in the beginning of the experimental curve, as the effect of dipolar fields vanishes towards the equilibration of the system. On Figure S11 the relaxation curves for DySc2N@C80 single crystal are shown. In the coordinates ued the square-root time dependence should appear as straight line, and one can see that only few first points on the experimental curve can be tentatively attributed to a square-root time dependence. Consequently, this specific regime in DySc2N@C80 is not accessible by the available magnetometry techniques.  Upon the increase of the relaxation time it becomes more and more time consuming to sample the complete relaxation curve. For the characteristic times of 10 5 seconds that would correspond to a 10 day experiment, which is already above the limit for a QD MPMS3 machine with unavoidable weekly liquid helium fillings (not mentioning the waste of the valuable SQUID time). The vital question then arises: how long is it necessary to measure the relaxation curve in order to obtain a reliable characteristic time?
The relaxation curve of DySc2N@C80 single crystal at 1.8 K, 70 mT is taken as an example. It has been measured for 2.8•10 4 seconds (~ 7 hours). Still in this case only a part of the decay curve was recorded ( Figure S12). It is especially problematic to determine the equilibrium magnetisation for such incomplete curves. In order to increase the reliability of the fit, an additional "thermal" magnetisation curve can be recorded. With the magnetic field kept at the initial value (70 mT in the current example), one warms the sample above the blocking temperature and then cools it down to the initial temperature, where the measurement starts. The magnetisation is rising towards the same Meq as for the decay curve. A combined fitting potentially yields more accurate results. In order to prove it, different fitting procedures are compared.
Single exponential fitting underestimates the characteristic relaxation time due to wrong value of Meq, as it can be seen on Figure S12. Once a second "rising" experimental curve is added to the single exponential fitting ("Single combined"), the obtained relaxation time becomes comparable to ones from the fittings with stretched exponentials. The shape of the simulated curve, however, is still far from experimental. This can be clearly seen at a closer look to the experimental curve ( Figure S13). Only stretched exponentials that can closely follow the curvature of the magnetisation decay curve with a time-dependent relaxation rate. For the full experimental curve the different stretched exponential fittings provide similar simulated curves and the relaxation times. Thus, the most reliable average relaxation time of (1.22±0.16)•10 5 seconds can estimated from 4 different fitting types: "Single combined", "Stretched", "Stretched combined", "Stretched combined normalized" (Figure S14). The ratio between the duration of the measurement and the characteristic time is ~0.2. Further, the cases with smaller ratios were modelled.
The initial curves (both decay and magnetisation) were continuously shortened from the end in steps of 2.8•10 3 seconds (10% of the full measurement). The resulting set of curves was fitted with a range of methods, both with default optimization algorithm and with a "normalised" one. However, the results for "normalised" fitting were excluded if they are similar to the default one ( Figure S14).
Following the Figure S14 one can see how the estimated relaxation time varies with the length on the experiment. In the case of single exponential fitting the estimated value changes continuously with the duration of the experiment. This is an additional argument against using this method.  The value of the relaxation time obtained by using "Single combined" fitting stays nearly constant down to ~0.1 ratios (experiment time to relaxation time), then rapidly deviates and occasionally coincides with single exponential fitting.
"Stretched" fitting values start to deviate immediately upon shortening of the experimental curve and the deviation can reach several orders of magnitude. This happens due to unconstrained Meq parameter that is difficult to estimate for a single decay curve.
"Stretched combined" method provides quite stable results, however, the standard deviation of the relaxation time value grows with shortening of the experiment. The reason is in the inappropriate weighting, which ignores the second "rising" magnetization curve as it has one order of magnitude smaller values than the decay curve ( Figure S12).
The weight of the data points from both decay and magnetisation experimental curves have to be normalized in the same way as it was done above for the case of short relaxation times. Accordingly, the "Stretched combined normalized" fitting method shows more accurate results, which do not deviate down to experiment time to relaxation time ratio of ~0.05 for the current case. Thus, the reliable experiment can't be more than 20 times shorter than the estimated relaxation time. The relaxation times determined for DySc2N@C80 samples in this work and listed in Tables S3-S7 were determined by this method.
Thus, we found that the stretched exponential is the optimal function for fitting the experimental DC relaxation curves. The following important notes can be pointed out:  Normalization of the weights is required when the measured signal varies within more than one order of magnitude.  For short relaxation times, the reliable values can be obtained down to ~100 s.  For long relaxation times, when only a part of the decay curve is measured, an additional measurement of the "rising" magnetization curve is required for an accurate estimation of Meq value. In this case two 30 h measurements will provide the reliable results for relaxation times up to ~10 6 seconds.   The function does not provide a good fit (deviations are significant). Besides, Raman process dominates at higher temperature and Orbach at lower temperature.  Figure S17. Relaxation times of magnetization of DySc2N@C80 (T>5 K, AC magnetometry), combination of zero-field measurements (T> 20 K) and in-field measurements (5 K < T < 20 K.) Fit with the combination of Orbach and Raman mechanisms similar to Fig. S16 above, but the exponent in the Raman process is also allowed to vary. As can be seen from the graph and from the errors in fitted parameters, the Orbach contribution in the fitted curve is negligible and hence non-reliable. Except for the two lowest-temperature points, the data are reasonably described by the power function of temperature with the exponent close to 3. This value is very different from the expected n = 9 for the Raman process in Kramers ion.  Figure S18. Relaxation times of magnetization of DySc2N@C80 (T>20 K, zero-field measurements, AC magnetometry). Fit with the combination of Orbach and Raman mechanisms, the exponent in the Raman process is also allowed to vary. The data is reasonably described by the power function of temperature with the exponent close to 2.4 (which is very different form the expected n = 9 for the Raman process in Kramers ion) and the Orbach process at higher temperatures. Note dramatic changes in the fitted parameters compared to the Fig. S17 caused by a reduction of the data set by few low-temperature points. This points to the low stability of the fit and hence makes the results questionable.

Ab initio calculations for DySc2N@C80
Single point CASSCF/RASSI calculation was done for a DFT-optimized structure of YSc2N@C80 (Table S9) with substitution of Y by Dy. For Dy 3+ ion the structure of high spin multiplet 6 H15/2 can be represented in eight low-lying Kramers doublets (Fig. S19). For ab initio modeling of the Dy multiplet structure we used Molcas 8.0 code. 15 The active space of the CASSCF calculations includes eleven active electrons and seven active orbitals (e.g. CAS (11,7)). All 21 sextet states and 108 quartets and only 100 doublets were included in the state-averaged CASSCF procedure and were further mixed by spin-orbit coupling in the RASSI procedure. VDZ quality atomic natural extended relativistic basis set (ANO-RCC) was employed. The single ion magnetic properties and ligand field parameters were calculated with use of SINGLE ANISO module (Table S10). The ab initio ligand field parameters were used to construct a model Zeeman Hamiltonian in |J,mJ>-basis (Clebsch-Gordan decomposition). Based on this Hamiltonian, the transition probabilities between different states were estimated using PHI program. 16

Simulation of magnetic hysteresis curves from spin dynamics
System settings. Let's consider the 2-level system with the energy gap Δ and coupling term (Ω ≠ 0): Such limited Hamiltonian is often used in resonance fluorescence studies. [17][18][19] Having supplied this system Hamiltonian into the GKSL equation, one will get a phenomenological master equation: By setting and the operator , we will allow one state to "relax" onto another with relaxation rate . We can set the environment coupling operator to be proportional to a raising (̂+) or lowering (̂−) operator depending on the relaxation direction. For example, if the lowering is chosen and ℒ(̂−) , the master equation would be read as following: Any generic two-level quantum system has the Hilbert space isomorphic to a spin (or pseudo spin)-half system and thus can be considered as one. The Hamiltonian Eq. S3 and the master equations Eq. S4/S5 together produce a suitable theoretical framework for magentodynamcis modeling. In that case the states energy split is linear with field Δ = Δ( ) ∼ and the parameter Ω is a non-zero probability to tunnel between these states. The expectation value of ⟨̂( )⟩ is magnetization by definition (⟨| ↑⟩⟨↑ |( ) − |↓⟩⟨↓ |( )⟩ ≡ ( )). Having that said, one can conduct a magnetization measurement experiment for observation period and in field H. To make this experiment possible, however, we have to introduce a self-evident field sweeping axioms.

Field sweep axioms
We assume that during the field sweep the system response is non-adiabatic locally, e.i. between two consequent field points, and adiabatic globally.
Let us prepare (or magnetize) the system Eq. S3 into the up-state by using field = ℎ. Next, we will sweep the field backwards from ℎ to -ℎ following the given rules:  from point to point the field is changed instantly, so is the Hamiltonian components (Δ( ) to Δ( +1 ));  at each point system stays for an observation period , during which magnetic moment evolved unperturbed according to equation of motion;  magnetic moment determination (observation) collapses the states on a "median-state", so that initial wavefunction operating on the new Hamiltonian (step: + 1) is defined as |Ψ⟩(0) = ⟨ ( )⟩ (0)| ↑⟩ + ⟨ ( )⟩ (0)| ↓⟩, where ( ) and ( ) expectation values for the operators | ↓⟩⟨↓ | and | ↑⟩⟨↑ | respectively, during observation time on the Hamiltonian at step ; Here is how this propagation algorithm works. Figure S20 shows for forward (in blue) and backward (in red) propagation of the systems described by Eq. S4. The dependence of time-averaged magnetization as a function of H or simply ℳ( ) ≡ ⟨ ( )⟩ is derived for different values of Ω in the range of 0.1-0.4 [ ℎ ], while the maximum energy difference is in the highest field Δ(ℎ) = 1[ ℎ ]. Thus, it appears that if Ω/Δ(ℎ) < < 0.1 on the field sweep from ℎ to −ℎ the state would not switch. Same field sweep algorithm is used in hysteresis modeling in the main text.

S33
Prove of correct thermodynamic limit and connection to LLG theory.
In the limit >> Ω for any and > 0 the master equation, so eq. S5 is read as: here ( ) is Boltzmann probability to find system in the up-state at temperature T.