Design of high-temperature f-block molecular nanomagnets through the control of vibration-induced spin relaxation

An efficient general first-principles methodology to simulate vibration-induced spin relaxation in f-block molecular nanomagnets that drastically reduces the computation time.


Introduction
Molecular nanomagnets have been attracting enormous attention for almost three decades due to their unique properties. These magnetic molecules, also known as single-molecule magnets (SMMs), 1 exhibit a bistable ground state that results in a memory effect characterized by a hysteresis loop and have the potential to be harnessed for classical information storage. The second generation of SMMs, 2 commonly called single-ion magnets (SIMs), is based on coordination complexes with a central magnetic ion as the source of magnetic anisotropy and represents the ultimate limit of miniaturization. Two key parameters are used to characterize the performance of SMMs and SIMs, namely, magnetic relaxation time and blocking temperature. The former is the timescale in which molecular nanomagnets preserve classical information at a given temperature, and the latter gives the maximum temperature that allows observing magnetic hysteresis.
The eld of molecular nanomagnetism is at a crucial point. Over the last few years, we have witnessed the discovery of new SIMs that have allowed an outstanding increase in the blocking temperature, rst from 30 K to 60 K (2017) 3,4 and then up to 80 K (2018). 5 This trend is drawing unprecedented interest from many researchers around the world and requires urgent attention from the theoretical point of view. Indeed, the large increase of the blocking temperatures opens the possibility to incorporate molecular nanomagnets in devices operating at higher temperatures. But rst, we need to deepen our understanding of the magnetic-behavior destroying process known as spin relaxation.
The target features that have commonly been addressed to block spin relaxation are (i) the ground electron spin quantum number J and (ii) the barrier height that separates the two components of the bistable ground state. 6 Indeed, the search for new nanomagnets by increasing (i) and (ii) is consistent with an Orbach relaxation mechanism, which drives the spin population across the barrier. While this strategy has worked up to now, 7,8 the recent interest in SIMs operating at high temperatures, where spin-vibration coupling dominates over relaxation, makes this scenario insufficient. 6 Hence, to gain control of relaxation at increasing temperatures, spinvibration coupling must be incorporated in the theoretical methods.
The current pursuit of predictive power is encouraging the development of fully ab initio methodologies. 9,10 Nevertheless, rst-principles evaluations of spin-vibration coupling are known to be computationally demanding. 3,9,10 This fact constitutes a crucial limitation that makes state-of-the-art ab initio methods impractical when guiding efforts at the lab stage. Thus, searching for new methodologies able to circumvent this computational bottleneck is of paramount urgency. In the case of lanthanide-based SIMs, there already exist affordable semiempirical approximations devoted to determining the electronic structure, 11,12 whose accuracy can become comparable to that of ab initio calculations. 13,14 Herein, we present an inexpensive rst-principles method devoted to lanthanide and uranium SIMs, with the aim of evaluating vibration-induced spin relaxation. It allows estimating effective relaxation barriers U eff , relaxation pathways and relaxation times s as a function of temperature. Crystal eld parameters (CFPs) are determined by millisecond calculations, and only one CASSCF evaluation is required. The method identies those vibrations promoting relaxation in order to redesign the given molecule and incorporates, for the rst time, a temperature dependence in the spin-vibration coupling matrix elements. Contributions from spin-spin dipole coupling to U eff and s can be incorporated by resorting to recent rstprinciples models. 15 Since the barrier height may increase from lanthanides to actinides due to a stronger ligand eld, and given the challenging computational nature of the U 3+ ion, 16,17 we propose to evaluate the effectiveness of bis-metallocenium ligands on actinides and test the efficiency of our method on the hypothetical analog [U(Cp ttt ) 2 ] + of [Dy(Cp ttt ) 2 ] + , which holds the latest record for the blocking temperature. 3 To assess the validity in SIMs with ligands of a different nature, we also apply this method to the experimentally studied uranium-based SIM UTp 3 (Tp À ¼ trispyrazolylborate) in order to rationalize its poor performance. The method consists of the following three steps.

Methods
Step 1 The relevant atom set is relaxed until reaching a minimum in its potential energy surface. 6 This set may be the SIM itself 9 or the unit cell of a crystal containing the SIM. 10 Aer calculating the vibrational spectrum, we extract harmonic frequencies {n j } j , reduced masses {m j } j , and displacement vectors fw j ! g j . These determine the 3D-direction in which each atom vibrates around its equilibrium position.
Step 2 We perform a CASSCF calculation on the SIM experimental structure to extract the lowest 2J + 1 energies, where J is the metal ion ground electron spin quantum number. Then, once the coordinate origin is placed at the metal ion experimental position, the crystallographic coordinates of only the metal-coordinating atoms are introduced in the code SIMPRE; see the ESI. † 11,12 This code rst calculates the CFPs by considering each coordinating atom to be an effective point charge and then performs a millisecond diagonalization of the ground J crystal eld Hamiltonian. We apply the Radial Effective Charge (REC) model by varying the magnitude of the effective charges and their radial distance to the metal ion to t the CASSCF energies; 18,19 see the ESI. † Thus, we project the CASSCF information onto the rst coordination sphere via effective REC parameters. Note that the contribution of the coordinating atoms to the ligand eld almost recovers the whole effect of magnetic anisotropy. Nevertheless, one can include non-coordinating ligand atoms if a signicant contribution is expected. To reduce computational costs, quite oen it will be enough to maintain the same charge magnitude Z i and radial distance contraction D r in each coordinating atom.
This procedure is fully ab initio, but one can avoid the CASSCF evaluation and use the experimental energies if they are available. In this case, the experimental structure used in SIMPRE should be determined at the same temperature as that of the experimental energies. Now, the coordinating atom positions of the relaxed geometry are radially varied with the same tting distance variations determined using SIMPRE. By using the same obtained charge values, SIMPRE calculates the equilibrium CFPs {(A k q hr k i) eq } k,q in Stevens notation. Unlike lanthanides, excited states beyond the ground J multiplet may also inuence the low-lying electronic structure of actinide SIMs. Thus, in the case of U 3+ , to determine the charge magnitude and the radial distance variation, the energy tting must be replaced by a tting of the SIMPRE CFPs to either the CASSCF or experimental CFPs. Yet, the energy tting procedure can be maintained providing SIMPRE is replaced with the package CONDON at the stage when the energy set is determined for each set of CFPs calculated using SIMPRE from the varying Z i and D r . CONDON contains the excited states beyond the ground J multiplet, and the CFPs must be introduced in Wybourne notation. 20 The diagonalization in SIMPRE of the equilibrium crystal the Stevens coefficient andÔ k q is the Stevens equivalent operator, 11,12 provides the lowest 2J + 1 equilibrium eigenstates and energies; see the ESI. † For U 3+ , the diagonalization is performed in CONDON. The lowest 2J + 1 equilibrium eigenstates obtained using this code are truncated to the ordered basis set {|ÀJi, ., |+Ji} of the ground J multiplet and then renormalized. The perturbed Hamiltonians , which are also built on the above ordered basis set, account for the perturbation to the equilibrium electronic structure from each vibrational mode j; see the ESI. † Their determination requires the estimation of the temperature-dependent change D(A k q hr k i) j (T) produced in (A k q hr k i) eq aer activating each mode j. We use a model derived by us elsewhere, 9 which provides the following perturbative expression up to the second-order in mode coordinate Q j : Thus, eachĤ j allows determining the spin-vibration coupling matrix element hi|Ĥ j |f i; see the ESI. † Temperature dependence is introduced for the rst time in these elements through each boson number hn j i ¼ 1/(e hn j /k B T À 1). The procedure to calculate the derivatives (v 2 A k q hr k i/vQ j 2 ) eq can be found in the ESI. † In any case, proper tuning of the chemical structure aimed at reducing the variations in D(A k q hr k i) j would improve the molecular nanomagnet performance. 9 Indeed, the aforementioned reduction would make the matrix elements hi|Ĥ j |f i and the transition rate g smaller, since the perturbed Hamiltonians Ĥ j are proportional to D(A k q hr k i) j . For instance, this could be undertaken by increasing the harmonic frequencies n j of the relevant vibrational modes with proper ligand modication. Of course, a legitimate question is whether these inexpensive calculations produce qualitatively similar CFPs at the equilibrium and distorted geometries as compared to the ones derived from ab initio calculations. This is what has been conrmed in a very recent study on the dysprosium-based SIM Dy-5*, 21 which holds the latest record for the blocking temperature.
Step 3 This step is undertaken by solving the master equation, eqn (2), 3,22,23 which describes the time evolution of the spin population across the lowest 2J + 1 equilibrium eigenstates. The energy that induces the spin to relax comes from the coupling with surrounding vibrations. Intuitively, at each time t there is a probability p i (t) of being in an eigenstate |ii. At a time t + dt the system may make a transition to a different eigenstate |f i with a probability g if dt, either by absorbing or by emitting a vibration quantum. The net difference between the incoming g  p i and outgoing g if p f spin populations equals the change in time of p i . Thus, once the transitions are assumed to be independent, these probabilities evolve as follows: The transition rates g if and g  account for the spin population ow between |ii and |fi, and their expressions depend on the relaxation process to model. We include two important processes: (i) Orbach and (ii) second-order Raman relaxation; see the ESI. † The determination of the most likely relaxation pathway provides further insight into relaxation. This allows identifying the vibrations that promote each relaxation step, and modications on the molecular structure can then be proposed to suppress relaxation. All details are found in the ESI. † 3,22-24 The same procedure applies if the crystal eld Hamiltonians are expanded to include a static magnetic eld via a Zeeman term.

Two case studies: UTp 3 and [U(Cp ttt ) 2 ] +
A handful of actinide complexes, mostly based on U 3+ , have been reported over the past few years as SIMs, although mostly with a poor magnetic behavior. 25 We choose to apply our theoretical approach to two complementary case studies. In the rst place, we studied UTp 3 , Tp À ¼ trispyrazolylborate, where the magnetic ion is directly bonded to nine pyrazole rings in a nearly exact D 3h tricapped trigonal prism coordination environment. 26 This is one of the few known uranium SIMs for which spectroscopic characterization has been performed and which thus can serve to validate our calculation of energy levels and wave-functions in these types of challenging systems. Secondly, we studied [U(Cp ttt ) 2 ] + , a hypothetical molecule that is analogous to the dysprosocenium SIM [Dy(Cp ttt ) 2 ] + . In this system, the f-ion is sandwiched between two tert-butyl(cyclopentadienyl) (Cp ttt ) ligands, and this gives rise to an overall linear coordination geometry that is slightly bent. Recently, a uranium-based SIM with a metallocenium-like ligand set similar to Cp ttt has been synthesized and experimentally characterized. This demonstrates the chemical viability of extending this particular family of ligands to uranium. 27 Since we lack the experimental structures of [U(Cp ttt ) 2 ] + and UTp 3 , we use the ones of [Dy(Cp ttt ) 2 ] + and NdTp 3 instead. 3,26 By replacing the corresponding lanthanide ion with U 3+ , we carry out geometry relaxation and vibrational spectrum calculation in both systems; see Fig. 1 and the ESI. † Now, we perform a CASSCF evaluation on real structures of [U(Cp ttt ) 2 ] + and UTp 3 to obtain {(A k q hr k i) eq } k,q . In the case of [U(Cp ttt ) 2 ] + , this is not possible and we will proceed as explained in the ESI † to determine the equilibrium electronic structure; see Fig. 1. Concerning UTp 3 , we will use the spectroscopic energy scheme to obtain its equilibrium electronic structure; see also the ESI † and Fig. 1.
In contrast with that of [U(Cp ttt ) 2 ] + , the energy level scheme of UTp 3 , see Fig. 1, is rather typical when compared with those of other studied uranium complexes. 17 The ground doublet of UTp 3 presents a heavy mixing between the |m J i components AE5/ 2 and H7/2, and this results in an almost perfect cancellation of the expectation value h J z i ¼ AE0.06. This seems to be a generalized feature that results from the combination of: (a) a nonperfect axial coordination geometry, where the AE5/2 and H7/2 components are favoured to the detriment of the maximum values AE9/2, and (b) the D 3h coordination symmetry, which results in a heavy mixing of levels differing by Dm J ¼ 6 (this is precisely the case for AE5/2 and H7/2). On the other hand, [U(Cp ttt ) 2 ] + also displays noticeable mixing as a consequence of the C 1 symmetry, although in this case the axial ligand distribution stabilizes a ground doublet mainly characterized by the component AE9/2 with a weight around 80%. This difference in the electronic structure already means that [U(Cp ttt ) 2 ] + should present a better prospect for slow relaxation of the magnetization as compared with UTp 3 .
Once (A k q hr k i) eq and (v 2 A k q hr k i/vQ j 2 ) eq are determined, we calculate the CFP thermal evolution. 9 From Fig. 2 Simple symmetry arguments can be used to give some intuitive meaning to these numerical results. In the case of [U(Cp ttt ) 2 ] + , the molecular symmetry is not ideal, meaning no CFPs cancel for symmetry reasons. As a consequence, there are no signicant changes caused to geometrical distortions of any reasonable size. This, added to the molecular rigidity that is characteristic of metallocenium complexes, results in a relative insensibility to thermal effects. Since the molecular geometry in this case is quite axial, the B 2 0 parameter remains dominant at all temperatures. The opposite situation happens for UTp 3 , where the perfect D 3h symmetry means many CFPs are cancelled at the equilibrium geometry. Thus, even moderate geometrical distortions cause dramatic relative changes in many CFPs, as indeed is the result of the calculations. In the special case of B 6 6 , this CFP is governed by a spherical harmonic of D 6h symmetry in such a way that all 9 donor atoms contribute with the same sign. Because the molecule presents an overall C 3 symmetry, certain concerted geometrical distortions that result from molecular vibrations couple with an unusually strong strength with this CFP: B 6 6 is at the same time the dominant parameter in the crystal eld Hamiltonian and the one with the largest absolute thermal effect. All in all, the dynamical effects coincide with the static picture in discarding the potential as a SIM of UTp 3 in particular, and possibly of C 3 -type U 3+ complexes in general.
In order to illustrate step 3 of the methodology we will now focus on the most promising case, namely [U(Cp ttt ) 2 ] + . Unfortunately, to the best of our knowledge, the uranium-based molecular nanomagnets reported so far in the literatureincluding UTp 3do not exhibit slow relaxation of the magnetization above ca. 5 K. This prevents us from applying step 3 to UTp 3 , since around and below this temperature a very high numerical noise leads to poorly reliable results; see the ESI. † The thermal dependence of the relaxation time s when the Orbach transition rates are used in eqn (2) is shown in Fig. 3. Above 30 K, where there exist a high enough number of available phonons, the thermally activated regime is at play and the spin population crosses the barrier through excited doublets; Fig. 4. Below 30 K, little or rather negligible spin population is promoted to the second excited doublet, which mostly tunnels the barrier through the rst excited doublet. Nevertheless, let us recall that as the temperature is decreased other mechanisms could come into play and even dominate over the Orbach-based one such as quantum tunneling between the ground doublet components. This process is not considered by our approach but commonly observed at low temperatures. 4,5 In this thermally activated regime, the estimated effective barrier U eff ¼ 293 K is in the range of hundreds of kelvin, which is normal in a large set of molecular nanomagnets, 28 and is  found around 40 cm À1 above the rst excited doublet in Fig. 1. The Orbach prefactor, s 0 ¼ e À19.62 ¼ 3.0 Â 10 À9 s, is within the usual range (10 À6 -10 À10 s) for SIMs with a barrier of a comparable height. We also evaluated eqn (2) with the second-order Raman transition rates. The Raman-based s values, see Table  S1, † are much larger than the ones in Fig. 3. Thus, the secondorder Raman process should be discarded as a competitive mechanism in the experiment, as found in [Dy(Cp ttt ) 2 ] + . 3 Even employing bis-metallocenium ligands, known to offer a strong axial crystal eld in Dy-based SIMs, 3-5 our calculated [U(Cp ttt ) 2 ] + effective barrier (293 K) is one order of magnitude below those reported for [Dy(Cp ttt ) 2 ] + (1760 K) 3 and Dy-5* (2217 K). 5 Besides, the maximum temperature in [U(Cp ttt ) 2 ] + ($40-50 K) at which the experimental relaxation time is still above the standard experimental detection limit (10 À5 -10 À6 s) is also clearly smaller compared to that in [Dy(Cp ttt ) 2 ] + ($112 K) 3 and Dy-5* ($138 K). 5 As can be expected for a decrease of about one order of magnitude in the barrier height, the calculated [U(Cp ttt ) 2 ] + Orbach prefactor s 0 ¼ 3.0 Â 10 À9 s is two to three orders of magnitude above the ones corresponding to the two Dy-based SIMs (s 0 $ 2.0 Â 10 À11 s and s 0 $ 4.2 Â 10 À12 s, resp.). 3,5 Our [U(Cp ttt ) 2 ] + Orbach prefactor is among the smallest ones that have been experimentally determined in uranium SIMs. 16,17 On the other hand, there do exist two signicant advances with respect to previous uranium SIMs: (i) the standard effective barrier is in the range of dozens of kelvin, 16,17 while that for [U(Cp ttt ) 2 ] + reaches several hundreds of kelvin ($293 K); (ii) by assuming that the thermally activated regime dominates between 30 K and 50 K in [U(Cp ttt ) 2 ] + , while it is not possible to determine relaxation times beyond $5 K, 16,17 the experimental detection limit s $ 10 À5 -10 À6 s in the case of [U(Cp ttt ) 2 ] + would be found at 40-50 K. Thus, it seems that [U(Cp ttt ) 2 ] + is not expected to display hysteresis temperatures that compete with the ones of [Dy(Cp ttt ) 2 ] + and Dy-5* SIMs. For instance, in a typical hysteresis loop swept at 2 mT s À1 between À1 T and +1 T, which means to store classical information for 2000 s, a blocking temperature below 9 K would be observed in [U(Cp ttt ) 2 ] + according to Fig. 3. This is unsurprising since, aer all, the Cp ttt rings were optimized for dysprosium and may present somewhat different requirements in uranium. However, our methodology is efficient enough to offer a path forward in the rational design of ligands that result in uranium SIMs with optimized performance.
Let us now analyze the [U(Cp ttt ) 2 ] + vibrations determining the transition rates that drive the relaxation pathways in Fig. 4, see the animations and the ESI, † and offer some strategies to reduce their detrimental effects. Two of them involve methyl rotations in the Cp ttt ring substituents. These rotations could be partially suppressed by replacing the methyl groups (-CH 3 ) with the heavier uorinated analogs -CF 3  has not yet been explored to remove the breathing vibrations could be bonding the two coordinating rings, such as in stapled bis-phthalocyanines. Moreover, the frequencies of the [U(Cp ttt ) 2 ] + detrimental vibrations, see the ESI, † closely match the gaps between the equilibrium ground and rst excited doublets (159.3 cm À1 ) and the rst and second excited doublets (171.7 cm À1 ); Fig. 1. Thus, the [U(Cp ttt ) 2 ] + performance would also benet from any structural modication that takes these frequencies out of resonance with respect to these gaps.

Conclusions
All in all, the most important conclusions of this work are the following. We have proposed a novel rst-principles methodology aimed to simulate vibration-induced spin relaxation in fblock SIMs. The method offers the important advantage of drastically reducing the computation time while keeping the calculation accuracy within an acceptable range. Indeed, all but one of the expensive CASSCF calculations required in previous methods are replaced by millisecond calculations. Besides, our approach introduces for the rst time a temperature dependence in the spin-vibration coupling matrix elements. To demonstrate this methodology, we consider two case studies. First, we revisit UTp 3 , a uranium SIM that has been studied spectroscopically and therefore allows us to apply our methodology with the highest-quality data. Here we nd a heavy mixing in the ground doublet which results in an almost perfect cancellation of the expectation value hJ z i. Furthermore, the intense thermal dependence of the CFPs evidences a strong coupling between molecular vibrations and the spin energy levels. This means that both from the static and from the dynamic point of view UTp 3 is expected to display poor SMM behavior. In the second place, we study the high-performing SIM [Dy(Cp ttt ) 2 ] + and nd that the replacement of Dy 3+ with U 3+ does not result in an enhanced performance. Yet, [U(Cp ttt ) 2 ] + does seem to outperform all previously reported uranium SIMs. One of the critical factors that promote spin relaxation in [U(Cp ttt ) 2 ] + is the noticeable mixing among the |m J i components in the equilibrium eigenstates. Importantly, this mixing is also found in previously reported uranium SIMs, 16,17,29 but not in the cutting-edge Dy-based SIMs [Dy(Cp ttt ) 2 ] + and Dy-5* even though the ligand coordination is not strictly axial. 3,5 Among those vibrations that promote spin relaxation, there are still atomic movements le to block. These involve methyl rotations and breathing deformations, which could be removed by uorination and stapling the coordinating rings to each other, respectively. Hence, we conclude that there may still be room for further improvement in these bis-metalloceniumbased uranium SIMs.

Conflicts of interest
There are no conicts to declare.