Deciphering the origin of giant magnetic anisotropy and fast quantum tunnelling in Rhenium(IV) single-molecule magnets

Single-molecule magnets represent a promising route to achieve potential applications such as high-density information storage and spintronics devices. Among others, 4d/5d elements such as Re(IV) ion are found to exhibit very large magnetic anisotropy, and inclusion of this ion-aggregated clusters yields several attractive molecular magnets. Here, using ab intio calculations, we unravel the source of giant magnetic anisotropy associated with the Re(IV) ions by studying a series of mononuclear Re(IV) six coordinate complexes. The low-lying doublet states are found to be responsible for large magnetic anisotropy and the sign of the axial zero-field splitting parameter (D) can be categorically predicted based on the position of the ligand coordination. Large transverse anisotropy along with large hyperfine interactions opens up multiple relaxation channels leading to a fast quantum tunnelling of the magnetization (QTM) process. Enhancing the Re-ligand covalency is found to significantly quench the QTM process.

Despite several years of comprehensive experimental efforts in designing Re(IV) ion-based SMMs/SCMs, the origin of giant magnetic anisotropy is not well understood 39 . As in most of the cases, the axial ZFS (D) values are extracted using magnetization measurements, which are known to be insensitive to the sign and strength of the D parameter 41,[43][44][45][46][47][48][49][50][51] . On other hand, the most promising high frequency-electron paramagnetic resonance (HF-EPR) technique has its own limitation, where the sign can be accurately determined but such large magnitude of D are often difficult to estimate 39 . An alternate solution to resolve the ambiguities in the sign/magnitude of D value is to analyse the ZFS parameter using ab initio calculations 9,22,39,[52][53][54][55][56][57] , which has been widely used in this respect. Moreover, strategic designing of new generation SMMs based on Re(IV) ions requires a thorough understanding of the nature ZFS and how the magnitude and the sign of the D and E vary depending on the ligand field environment. As the magnitude of E and the hyperfine interactions are correlated to the quantum tunnelling of magnetization (QTM) 58 , the possibility to fine tuning these values is of paramount importance in this area.
The goal of the present communication is to gain a thorough understanding of the magnetic anisotropy in six coordinate Re(IV) complexes using state-of-the-art ab initio calculations. By modelling structurally diverse 13 mononuclear six coordinate Re(IV) complexes 39,41,[43][44][45][46][47][48][49][50][51] , we aim to answer the following intriguing questions (i) What is the suitable theoretical methodology to compute ZFS parameters in 5d transition metal ions such as Re(IV) complexes? (ii) What is the origin of giant D values and is there a correlation between the nature of the donor atoms and the sign of the D values in these complexes? (iii) What is mechanism of magnetic relaxation in Re(IV) singleion magnets and how this is influenced by the metal-ligand covalency?

Results
Magnetic anisotropy and spin-Hamiltonian. The free Re(IV) ion is a d 3 Kramers ion with a 4 F ground state term, which splits into three states 4 A 2g , 4 T 2g and 4 T 1g with the 4 A 2g being the ground state in an octahedral environment. Due to perfect cubic symmetry, pure octahedral complexes do not possesses any ZFS, however, any distortions from the octahedral geometry are expected to yield large D values via mixing of the subsequent excited states because of very large SOC (lB1,000 cm À 1 ). To begin with, we have studied the homoleptic [ReCl 6 ] 2 À model complex in tetragonal environment to analyse the origin of ZFS. Abragam and Bleaney have proposed a qualitative equation to predict the sign of the D values of a tetragonally distorted d 3 ion.
where l is the SOC and D 0 , D 1 is the tetragonal distortion parameter. If D 0 oD 1 , the sign is predicted to be negative, whereas for D 0 4D 1 , the sign is predicted to be positive. To corroborate this qualitative analysis, we have performed ab initio calculations, using complete active self consistent field (CASSCF) and CASPT2 methods incorporating spin-orbit effects with the RASSI-SO module in MOLCAS 59,60   Sign and magnitude of ZFS parameter for 1-13. Calculations reveal that eight spin-free states corresponding to 2 G states are found to be low-lying and thus are expected to contribute significantly to the D values via spin-flip excitations in all complexes 1-13 studied (see Supplementary Figs 3   witnessed. On other hand, complexes of type II and type III categories (7)(8)(9)(10)(11)(12)(13) found to posses positive ZFS parameter with D as high as þ 55 cm À 1 has been noticed (for 8; see Supplementary  Tables 17 and  18 for details). For complex 6, the magnetization data 41 yield an estimate of D as À 14.4 cm À 1 , which is in agreement with CASSCF results ( À 21 cm À 1 , see Supplementary Table 17) but in disagreement with MS-CASPT2 values ( þ 16.2 cm À 1 , see Table 1). However, HF-EPR experiments performed lately 35 , where both the magnitude as well as the sign of the D value is estimated accurately, places the D value to be þ 11 cm À 1 . This highlights the issue of obtaining the sign/magnitude of the D value from the magnetization data and also emphasise the need for CASPT2 approach and hence incorporation of dynamic correlation to correctly reproduce the sign of the D values 62,63 . Inclusion of dynamic correlation on CASSCF computed wavefunctions drastically stabilizes the doublet states compared with the computed CASSCF states, leading to a pronounced contribution from these states to the D values as discussed earlier (see Supplementary Fig. 3 for details). Moreover, either pseudospin or effective Hamiltonian approach 62 needs to be employed to extract the ZFS parameters as other theoretical methodologies found to yield ambiguous sign and magnitude of D values (see Supplementary Tables 19 and 20 for details).
The SINGLE_ANISO computed orientations of the main anisotropic axes (D XX , D YY and D ZZ ) and main magnetic axes (g XX , g YY and g ZZ ) for all the complexes 1-13 are provided in   Fig. 3 and Supplementary Fig. 5. It is evident from the figures that, the principal anisotropic axes (D ZZ ) are oriented towards the L ax -R-L ax (molecular -z axis) direction, however, a significant tilt from this axis is witnessed 64 . Moreover, the main magnetic axes (g XX , g YY and g ZZ ) and main anisotropic axes (D XX , D YY and D ZZ ) do not coincide with each other and such non-coincidence has been previously noticed by Askevold et al. 65 The orientation of the D ZZ axis is tilted by 28 Fig. 4 for details). For 1, presence of unsymmetrical ligand in the equatorial position leads to the d xz orbital being the lowest lying in energy followed by degenerate d yz and d xy orbitals. The p* orbitals of the oxygen donors interact rather strongly with the d yz and d xy orbitals because of acute +O-Re-O bite angle. The d xz orbital on the other hand faces less repulsion leading to a slight stabilization. Besides, stronger s* interaction by the oxalate ligand lead to destabilization of the d x 2 -y 2 orbital compared with the d z 2 orbital. This orbital ordering has the following consequences to the D values: (i) all spinconserved excitations from the d xy , d xz and d yz orbitals to the vacant d z 2 orbital contribute to positive D values. This is affirmed by an additional calculation incorporating only the quartet states in the estimation of D and this yield a positive D value ( þ 13 cm À 1 ). (ii) Spin-flip excitations from the d xz to the d yz orbitals contribute to negative D values (excitation between same |m l | levels). As the gap between these two orbitals is very small (667 cm À 1 ), this transition governs both the sign and the magnitude of the D value for this complex. A similar pattern predicted for complexes 2-6 rationalize the observed negative sign for these complexes. Strong p acceptor cyanide ligand stabilizes the d xz and d yz orbitals via pp-dp interactions compared with the d xy orbital. Stronger s-donation along the axial directions destabilizes the d z 2   orbital leading to different orbital arrangement compared with complex 1. This orbital arrangement has the following consequences to the D value: (i) the spin-conserved d xy -d x 2 -y 2 transition contributes to the negative D value, affirmed by an additional calculation incorporating only the quartets in the estimation of D and this yields À 26 cm À 1 as the D value, (ii) the spin-flip excitations from the degenerate d xz -d yz orbitals to the d xy orbital contribute to the positive D value (the gap here is 3,123 cm À 1 ). Here also the later term dominates leading to an overall positive D value for 7.
For complex 10, positive D value is expected as the splitting pattern is found to be very similar to that of complex 7 (Fig. 4). However, the presence of weak p-donor pyridazine ligand on the axial position reduces the splitting between the d xz -d yz orbital and the d xy orbital (1,044 cm À 1 ). This suggests that much less energy is required to flip the spin in complex 10 compared with complex 7, thus the D value is expected to be large in this case. A similar splitting pattern is predicted for complexes 8, 9 and 11-13 rationalize the observed positive D values for these complexes.
Rationale for the observed variation in the ZFS parameter. Re-Br distance same as that of complex 2. For this model, the D is estimated to be À 86 cm À 1 compared with À 85 cm À 1 for the Cl analogue and this suggest that apart from the spin-orbit coupling of Br, the structural distortion such as -cis angles play an important role in determining the strength of the D value (see Supplementary Table 3 for selected structural parameters of complexes 1 and 2) 55 .
The strength of the donor-acceptor abilities significantly affects the magnitude of the D values (see Fig. 5 and Supplementary Table 21 for further details) 20,62 . In complexes 1-6, larger charge on the donor atoms are found to yield large D values (see Fig. 4 and Supplementary Figs 11 and 12 for quantitative charges computed). Larger charges on the donor atoms stabilizes the d xz orbitals compared with the d xy /d yz orbitals leading to different transition energies (see Table 1 for details) and thus the computed charges are found to strongly correlated to the magnitude of the D values. This striking observation offers a rational approach to fine tune the magnitude of the negative D value in this set of complexes. Moreover, stabilization of d xz orbital also affects the E values, as it increases the difference between the D XX and D YY contributions leading to a larger E with large charge on the ligand (see equation 2 and Supplementary  Fig. 12 for details).
In contrast to complex 7, where axial positions are occupied by two strong p-acceptor ligands, complex 8 posses two weak p-donor (pyridine) ligands on the axial position and therefore the d xz -d yz orbital to d xy orbital is found to be (1,360.7 cm À 1 ) much smaller than that of 3,123 cm À 1 observed in case of complex 7. Moreover, the first spin-flip-excited state is found at 5,950.5 cm À 1 (in case of complex 8), which is again much smaller than the 7,377.3 cm À 1 gap observed for complex 7 (see Table 1 for details). This leads to larger D value for complex 8 compared with complex 7. Complexes 9-13 possess positive D values ranging from þ 18 cm À 1 (complex 13) to þ 41 cm À 1 (complex 10). As the structural parameters across the series are very similar, the differences in the magnitude of the D values are expected to arise from the donor strength of the ligand. To affirm this point, we have analysed the donor-acceptor interactions using secondorder perturbation theory natural bonding orbitals (NBO) analysis for complexes 10 and 13. NBO analysis suggests a significant s-donation from lone pair of nitrogen to Re d z 2 orbital and this strength is estimated to be 19.6 kcal mol À 1 for complex 10, whereas 22.6 kcal mol À 1 for complex 13 (see Supplementary  Fig. 13 for details). Larger s-donation in complex 13 leads to smaller D value compared with complex 10. A similar analogy can be drawn also for other complexes.
Here independent of the nature of p-donor/acceptor ligands, type I complexes found to yield negative D values, whereas type II and III complexes found to yield positive D values. This is in stark contrast to the earlier observations where lighter transition metal d 3 ions found to switch the sign of ZFS parameter by changing the nature of the ligand donor atoms (p-donor ligands found to yield þ D values, whereas p-acceptor ligands yield -D values) 66 . This is essentially due to the fact that spin-flip doublet transitions are the dominating factor to the D values in Re(IV) complexes, whereas in lighter elements due to smaller crystal-field splitting, the spin-allowed transition dominates the D value.
Magneto-structural D-correlations. To probe how the axial bond length influences the D value, we have developed a magneto-structural D correlations on complex 7 (see Fig. 5c  , where the axial -CN bonds are varied from 2.0 to 2.6 Å (compression and elongation of axial bonds). As [ReCl 4 (CN) 2 ] 2 À unit has been employed as a building block for synthesis of polynuclear SMMs/SCMs [34][35][36]41 , magneto-structural correlation developed on this model will serve the purpose of obtaining qualitative single-ion Re(IV) anisotropy in diverse polynuclear framework. This Re-C bond distance is found to vary significantly among structures, particularly when the -CN ligands are found to bind to other metal ions. In our correlation, the magnitude of the D value found to drastically increase (from þ 13.6 to þ 94.7 cm À 1 ) as the Re-C bond length increases to 2.7 Å. As the metal-ligand interactions are weaker at longer distances, the transition energies are further lowered leading to larger D values for axially elongated structures 62,63 .
Besides our results reveal that tetragonal distortions does not alter the sign of D values and this suggests that the [ReCl 4 (CN) 2 ] 2 À unit unlikely to offer negative single-ion D value in any polynuclear framework.
Mechanism of magnetic relaxation. Complex 1 exhibits fieldinduced SMM behaviour with a barrier height of 9.6 cm À 1 at higher temperatures and 1.5 cm À 1 at lower temperatures. The experimental relaxation observed at higher temperature (up to 3.5 K) is unlikely due to Orbach process as the first excited Kramer's doublet (KD) is estimated to lie at 195 cm À 1 . Thus, the relaxation is expected to be a multi-phonon Raman process. The fast relaxation observed at lower temperature is essentially due to QTM process, which is facilitated due to the presence of transverse anisotropy, hyperfine interactions and external perturbations such as internal magnetic field provided by surrounding molecules. To gain insights into the QTM process, we have analysed the wavefunction of the ground-state KD and our analysis suggests that the ground-state KD comprised of 44% of |3/ 2, ± 3/2i and 47% of |3/2, ± 1/2i. As the D value is very large (D44kT), the |3/2, ± 1/2i KD will be completely depopulated and the ground state can be treated as a pseudo spin 1/2 system. The presence of large E term offers a strong mixing between the |3/2,±3/2i and |3/2,±1/2i components, which allows QTM to facilitate at low temperatures. To qualitatively analyse the mechanism of magnetic relaxation, we have computed the matrix elements between the connected KDs (see Supplementary Fig. 14 and Supplementary Note 3 for details) 67 . Our calculations predict very large tunnelling probability between the ground-state KDs and this is in line with the analysed wave function analysis. Such a prominent QTM process expected to quench the magnetization completely and this is consistent with the absence of zero-field SMM behaviour 56

Discussion
Here we have probed the origin of contrasting behaviour observed for Re(IV) SMMs where both the giant magnetic anisotropy and fast QTM found to co-exist. The low-lying doublet states are found to govern the sign and magnitude of ZFS parameters in this class of complexes. Our method assessment reveals that pseudo spin approach or effective Hamiltonian approach coupled with CASPT2 calculations needs to be employed to correctly reproduce the sign and magnitude of ZFS parameters. Quite interestingly, the sign of D values are found to be predictable based on the coordination mode of the ligands in these complexes where type I complexes found to possess larger negative D values, whereas type II and III found to possess a positive D values. Nature of the donor ligands as well as charge on the coordinated atoms found to influence only the strength but not the sign of D values. Very large hyperfine interactions (both transverse and axial) and rhombic anisotropy computed on these system found to govern the QTM process. By performing additional calculations and by developing magnetostructural correlations, we offer a way to enhance (diminish) the negative D (|E|) value in these classes of complexes. Some interesting observations are noted where the metal-ligand covalency found to govern the transverse anisotropy, offering a way to quench the inherent fast QTM process in this class of complexes.

Methods
Ab initio calculations. We have performed the ab initio calculations based of wave function theory approach to compute the ZFS in these set of mononuclear complexes. All the calculations have been performed using MOLCAS 7.8 suite of programme 69 . Here we have employed the state average-CASSCF method to compute the ZFS. The active space comprises of three active electrons in five active orbitals (CAS (3,5)). With this active space, we have computed all the 10 quartets and 40 doublet states in the configuration interaction procedure. On top of the converged CASSCF wave function, we have performed MSCASPT2 calculations to treat the dynamical correlations. We have employed ionization potential electron affinity (IPEA) shift of 0.25 to avoid the intruder states problem in CASPT2 calculations. The MS-CASPT2 computed states were further treated in RASSI-SO module, which explicitly computes the spin-orbit states. Furthermore, SINGLE_ANISO module has been utilized on top to compute the reliable spin-Hamiltonian (D and g-tensor, orientation of main magnetic axes and main anisotropic axes and local magnetic susceptibility) for each complex. DFT calculations. Hyperfine interaction of the Re(IV) nuclei were computed within DFT framework, using electron paramagnetic resonance/nuclear magnetic resonance (EPR/NMR) module in the ORCA code 70 . We have employed meta-GGA TPSSH functional along with SARC basis set for the Re, which is much more flexible at core region to estimate all the components of the A-tensors (Fermi Fermi Contact, Spin-dipolar and Spin-orbit coupling) along -x, -y and -z directions. A very tight self consistent field (SCF) (1 Â 10 À 8 E h ) has been kept throughout the calculations.