Spin-lattice relaxation of magnetic centers in molecular crystals at low temperature

TL;DR

This paper derives universal, measurable formulae for spin-phonon relaxation rates in molecular magnets at low temperature, revealing dependencies on magnetic field orientation and temperature, and proposes strategies for improving single-molecule magnet performance.

Contribution

The paper introduces universal, measurable formulae for spin-phonon relaxation rates that are independent of excited state energy gaps, applicable to both Kramers and non-Kramers systems.

Findings

Relaxation rates expressed via measurable quantities

Raman process depends on magnetic field orientation

Raman process varies as T^9 in both systems

Abstract

We study the spin-phonon relaxation rate of both Kramers and non-Kramers molecular magnets in strongly diluted samples at low temperature. Using the "rotational" contribution to the spin-phonon Hamiltonian, universal formulae for the relaxation rate are obtained. Intriguingly, these formulae are all entirely expressed via measurable or \emph{ab initio} computable physical quantities. Moreover, they are also independent of the energy gaps to excited states involved in the relaxation process. These obtained expressions for direct and Raman processes offer an easy way to determine the lowest limit of the spin-phonon relaxation of any spin system based on magnetic properties of the ground doublet only. In addition, some intriguing properties of Raman process are also found. Particularly, Raman process in Kramers system is found dependent on the magnetic field's orientation but independent of its magnitude, meanwhile the same process in non-Kramers system is significantly reduced out of resonance, i.e. for an applied external field. Interestingly, Raman process is demonstrated to vary as $T^{9}$ for both systems. Application of the theory to a recently investigated cobalt(II) complex shows that it can provide a reasonably good description for the relaxation. Based on these findings, a strategy in developing efficient single-molecule magnets by enhancing the mechanical rigidity of the molecular unit is proposed.

Equations

$$H_{ZFS}=D{S}_z^2+E(S_x^2-S_y^2),$$

$$H_{\text{CF}}=\sum _{p=2,4,6}\sum _{q=-p}^p{B}_{pq}{O}_p^q(\mathbf{J}),$$

$$\frac{d{\rho }_{\alpha \beta }^{\left(I\right)}}{dt}=\sum _{{\alpha }^{\prime },{\beta }^{\prime }}{e}^{i\left({\omega }_{\alpha \beta }-{\omega }_{{\alpha }^{\prime }{\beta }^{\prime }}\right)t}{R}_{\alpha \beta ,{\alpha }^{\prime }{\beta }^{\prime }}{\rho }_{{\alpha }^{\prime }{\beta }^{\prime }}^{\left(I\right)},$$

$${\mathcal{H}}_{\mathrm{total}}=R{\mathcal{H}}_A{R}^{-1}+{\mathcal{H}}_Z+{\mathcal{H}}_{\mathrm{ph}},$$

$${\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}$$

$${\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}^{\left(1\right)}$$

$${\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}^{\left(2\right)}$$

$$\delta \phi =\frac{1}{2}\nabla \times \mathrm{u}\left(\mathrm{r}\right)=\sqrt{\frac{\hslash }{8\rho V}}\sum _{\mathrm{k}{\lambda }_{\mathrm{k}}}\frac{\left[i\mathrm{k}\times {\mathrm{e}}_{\mathrm{k}{\lambda }_{\mathrm{k}}}\right]e^{i\mathrm{k}\cdot \mathrm{r}}}{\sqrt{{\omega }_{\mathrm{k}{\lambda }_{\mathrm{k}}}}}\left(a_{\mathrm{k}{\lambda }_{\mathrm{k}}}+a_{-\mathrm{k}{\lambda }_{\mathrm{k}}}^{+}\right).$$

$${\Gamma }_{-+}^{\mathrm{dr}}$$

$${\Gamma }_{-+}^{\mathrm{Raman}}$$

$$M^{\left(1\right)}\equiv ⟨-;n_{\mathrm{k}{\lambda }_{\mathrm{k}}}+1⟩{\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}^{\left(1\right)}+;n_{\mathrm{k}{\lambda }_{\mathrm{k}}},$$

$$M^{\left(2\right)}\equiv ⟨-;n_{\mathrm{k}{\lambda }_{\mathrm{k}}}-1,n_{\mathrm{q}{\lambda }_{\mathrm{q}}}+1⟩{\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}^{\left(2\right)}+;n_{\mathrm{k}{\lambda }_{\mathrm{k}}},n_{\mathrm{q}{\lambda }_{\mathrm{q}}},$$

$$M^{\left(1+1\right)}\equiv \sum _{\xi }\sum _{\begin{array}{c}{\mathrm{k}}^{\prime }{\lambda }_{{\mathrm{k}}^{\prime }}\\ {\mathrm{q}}^{\prime }{\lambda }_{{\mathrm{q}}^{\prime }}\end{array}}\sum _{\begin{array}{c}{n}_{{\mathrm{k}}^{\prime }{\lambda }_{{\mathrm{k}}^{\prime }}}\\ n_{{\mathrm{q}}^{\prime }{\lambda }_{{\mathrm{q}}^{\prime }}}\end{array}}\frac{⟨-;n_{\mathrm{k}{\lambda }_{\mathrm{k}}}-1,n_{\mathrm{q}{\lambda }_{\mathrm{q}}}+1⟩{\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}^{\left(1\right)}\xi ;n_{{\mathrm{k}}^{\prime }{\lambda }_{{\mathrm{k}}^{\prime }}},n_{{\mathrm{q}}^{\prime }{\lambda }_{{\mathrm{q}}^{\prime }}}⟨\xi ;n_{{\mathrm{k}}^{\prime }{\lambda }_{{\mathrm{k}}^{\prime }}},n_{{\mathrm{q}}^{\prime }{\lambda }_{{\mathrm{q}}^{\prime }}}⟩{\mathcal{H}}_{\mathrm{sp}-\mathrm{ph}}^{\left(1\right)}+;n_{\mathrm{k}{\lambda }_{\mathrm{k}}},n_{\mathrm{q}{\lambda }_{\mathrm{q}}}}{{\epsilon }_{+}-{\epsilon }_{\xi }+E_{n_{\mathrm{k}{\lambda }_{\mathrm{k}}},n_{\mathrm{q}{\lambda }_{\mathrm{q}}},\dots }-E_{n_{{\mathrm{k}}^{\prime }{\lambda }_{{\mathrm{k}}^{\prime }}},n_{{\mathrm{q}}^{\prime }{\lambda }_{{\mathrm{q}}^{\prime }}},\dots }}.$$

$$M^{\left(1\right)}=\sum _{\alpha }{\Xi }_{\alpha }{M}_{\mathrm{ph}}^{\alpha },$$

$${\Xi }_{\alpha }\equiv ⟨-⟩{\mathcal{H}}^{\left(1\right)\alpha }+,$$

$$M_{\mathrm{ph}}^{\alpha }\equiv \sqrt{\frac{\hslash }{8\rho V}}\frac{{\left[-i\mathrm{k}\times {\mathrm{e}}_{\mathrm{k}{\lambda }_{\mathrm{k}}}\right]}_{\alpha }{e}^{-i\mathrm{k}\cdot \mathrm{r}}}{\sqrt{{\omega }_{\mathrm{k}{\lambda }_{\mathrm{k}}}}}\sqrt{n_{\mathrm{k}{\lambda }_{\mathrm{k}}}+1}.$$

$${\Gamma }_{-+}^{\mathrm{dr}}$$

$$M^{\left(2\right)}$$

$$=\sum _{\alpha ,\beta }⟨-⟩{\mathcal{H}}^{\left(2\right)\alpha \beta }+{\mathcal{H}}^{\left(2\right)\beta \alpha }+M_{\mathrm{ph}}^{\alpha \beta },$$

$$M_{\mathrm{ph}}^{\alpha \beta }$$

$$=\frac{\hslash }{8\rho V}\frac{{\left(\mathrm{q}\times {\mathrm{e}}_{\mathrm{q}{\lambda }_{\mathrm{q}}}\right)}_{\alpha }{\left(\mathrm{k}\times {\mathrm{e}}_{\mathrm{k}{\lambda }_{\mathrm{k}}}\right)}_{\beta }}{\sqrt{{\omega }_{\mathrm{q}{\lambda }_{\mathrm{q}}}{\omega }_{\mathrm{k}{\lambda }_{\mathrm{k}}}}}\sqrt{\left(n_{\mathrm{q}{\lambda }_{\mathrm{q}}}+1\right)n_{\mathrm{k}{\lambda }_{\mathrm{k}}}}{e}^{i\left(\mathrm{k}-\mathrm{q}\right)\mathrm{r}}.$$

$$M^{\left(1+1\right)}$$

$$M^{\left(2\right)}+M^{\left(1+1\right)}$$

$$+\sum _{x_m}\left(\frac{{\mathcal{H}}_{m_1{x}_m}^{\left(1\right)\alpha }{\mathcal{H}}_{x_m{p}_1}^{\left(1\right)\beta }}{\hslash {\omega }_{+x_m}+\hslash {\omega }_{\mathrm{k}{\lambda }_{\mathrm{k}}}}+\frac{{\mathcal{H}}_{m_1{x}_m}^{\left(1\right)\beta }{\mathcal{H}}_{x_m{p}_1}^{\left(1\right)\alpha }}{\hslash {\omega }_{+x_m}-\hslash {\omega }_{\mathrm{q}{\lambda }_{\mathrm{q}}}}\right)]M_{\mathrm{ph}}^{\alpha \beta }\equiv \sum _{\alpha ,\beta }{M}_{\mathrm{sp}}^{\alpha \beta }{M}_{\mathrm{ph}}^{\alpha \beta },$$

$$\Omega \ll \hslash {\omega }_{\mathrm{k}{\lambda }_{\mathrm{k}}}\ll \hslash {\omega }_{\xi 1},|\xi ⟩\ne |\pm ⟩.$$

$${\mathcal{H}}_A=\sum _m{\epsilon }_m\left(|m^{\left(0\right)}⟩⟨m^{\left(0\right)}|+|{\bar{m}}^{\left(0\right)}⟩⟨{\bar{m}}^{\left(0\right)}|\right),$$

$${\left({\mathcal{H}}_A+{\mathcal{H}}_Z\right)}_{\left\{|1^{\left(0\right)}⟩,|{\bar{1}}^{\left(0\right)}⟩\right\}}$$

$${\Delta }_x=g_x{H}_x,{\Delta }_y=g_y{H}_y,$$

$$\Delta =\sqrt{{\Delta }_x^2+{\Delta }_y^2},W=g_z{H}_z,$$

$$\sin \phi =\frac{{\Delta }_y}{\Delta },\cos \phi =\frac{{\Delta }_x}{\Delta }.$$

$$|1⟩$$

$$|M⟩$$

$${\left({\mathcal{H}}_A+{\mathcal{H}}_Z\right)}_{\left\{|1⟩,|\bar{1}⟩\right\}}=\frac{W}{2}\left(|1⟩⟨1|-|\bar{1}⟩⟨\bar{1}|\right)+\frac{\Delta }{2}\left(|1⟩⟨\bar{1}|+|\bar{1}⟩⟨1|\right).$$

$${\epsilon }_{\pm }=\pm \frac{1}{2}\sqrt{W^2+{\Delta }^2}\equiv \pm \frac{1}{2}\Omega ,$$

$$|p_1⟩=\cos \frac{\vartheta }{2}|1⟩+\sin \frac{\vartheta }{2}|\bar{1}⟩,$$

$$|m_1⟩=\cos \frac{\vartheta }{2}|\bar{1}⟩-\sin \frac{\vartheta }{2}|1⟩=|{\bar{p}}_1⟩.$$

Full text

Single-molecule magnets, spin-phonon relaxation.

††preprint: 0APS/123-QED

I Introduction

The spin-lattice interaction plays a major role in a number of application-potential physical phenomena such as magnetic relaxation in single-molecule/single-atom magnets Gatteschi et al. (2006) (SMMs), entanglement of quantum system to the environment Horodecki et al. (2009), and decoherence of magnetic qubits Zurek (2003). In particular, for technological applications of single-molecule magnets, a low spin-lattice relaxation rate is a critical factorGatteschi and Sessoli (2003); Gatteschi et al. (2006) and a complete understanding of mechanisms governing this rate is a must. However, since the accurate evaluation of this spin-phonon relaxation rate requires knowledge of corresponding matrix elements of the deformation potential with respect to all phonons in the crystal Van Vleck (1940), a universal form of this rate is hardly achievable.

In the case of spin complexes, such as Mn12ac and Fe8 Gatteschi and Sessoli (2003), the magnetic relaxation proceeds via zero-field-split (ZFS) components of the ground state spin $S$. Due to quenched orbital momentum in the ground state of these complexes the energy separations between the ZFS components of the $S$-term do not exceed 10$\div$20 cm*-1*, i.e., they are in the range of frequencies of acoustic phonons even in molecular crystals. The latter are well described by the Debye model resulting in expressions for spin-phonon relaxation rates Abragam and Bleaney (1970) which have been widely employed for the description of magnetization relaxation in $S$-complexes Villain et al. (1994); Hartmann-Boutron et al. (1996); Garanin and Chudnovsky (1997); Leuenberger and Loss (2000). Besides the limitation to the acoustic phonons, the description of spin-phonon interaction requires knowledge of derivatives with respect to deformations of two parameters defining the ZFS of a spin manifold 111Besides modifying the ZFS parameters, the deformation also induces the rotation of the main anisotropy axes of the ZFS tensors described by three parameters (Euler angles) Chibotaru (2013) :

[TABLE]

where $S_{\alpha}$, $\alpha=x,y,z$, are spin operators acting on the spin states of the ground $S$ term. Often only the contribution of the parameter of axial anisotropy ($D$) is considered while the contribution of the parameter of rhombic anisotropy ($E$) is usually neglected due to its smallness in SMMs Garanin and Chudnovsky (1997); Leuenberger and Loss (2000). On the same reason are neglected also the quartic terms in the ZFS Hamiltonian (1) and their contribution to the spin-phonon interaction.

Despite this great simplification, the description of spin-phonon relaxation in SMM $S$-complexes still requires the knowledge of derivative of the parameter $D$ and $E$ and directions of the anisotropy axes Chibotaru (2013) with respect to all phonons. In connection with this problem, Chudnovsky et. al. Chudnovsky et al. (2005) made an important observation that given the strong dependence of the transition rate on the velocity of sound of acoustic phonons, the most efficient should be the relaxation via the transversal phonon modes. In molecular crystals, the deformations related to the transversal phonons and concomitantly strongly affecting the ZFS interaction (1) are the vibration modes corresponding to small rotational displacements of the undeformed SMM molecule from the equilibrium point. The variation of $H_{\mathrm{ZFS}}$ under these distortions reduces to the rotation of its main anisotropy axes following the rotation of the molecule without modifying the parameters $D$ and $E$. Given the kinematic nature of this “rotational” contribution to the spin-phonon interaction, the corresponding expression for the transition rate between two tunneling states is found to be solely defined by the form of the spin Hamiltonian and independent of any spin-phonon phenomenological coupling parameters Chudnovsky et al. (2005). This remarkable finding allows to evaluate the main contribution to the relaxation rate in $S$-complexes at low temperature without knowledge of the details of spin-phonon interaction. In a general case when other deformations play important role, a calculation of this rotational contribution to the relaxation rate is also meaningful in that it gives the lower boundary for the relaxation rate, which is a crucial insight for the design of efficient SMMs.

In lanthanide-based complexes, which is attracting increasing attention in the last years due to their exceptional magnetization blocking Layfield and Murugesu (2015); Woodruff et al. (2013); Ishikawa et al. (2003), the spin-phonon interaction is by far more complex. The reason is that the operator describing the crystal-field (CF) splitting of the ground $J$ atomic-multiplet of the lanthanide ion,

[TABLE]

where $O_{p}^{q}\left(\bm{\mathrm{J}}\right)$ are the extended Stevens operators Rudowicz (1985); Rudowicz and Chung (2004); Gatteschi et al. (2006), contains generally 27 parameters instead of five ($D$, $E$, and the directions of anisotropy axes) in the case of $S$-complexes. Taking into account the variations with local distortions of all parameters $B_{pq}$ results in a highly complex form of spin-phonon interaction Dohm and Fulde (1975), containing a tremendous number of parameters. On the other hand the arguments given by Chudnovsky et. al. Chudnovsky et al. (2005) should be applicable for lanthanide complexes too, for which the derivation of a simple universal expression for spin-phonon transition rate is especially desirable. Additionally, in lanthanides the ground CF doublet is usually well separated from the excited ones. Accordingly, the magnetic relaxation in such complexes takes place preponderantly within the ground doublet in a broad temperature domain and the relaxation reduces to the spin-phonon transition between its two tunnel-split states Garanin (2011).

Given the general complexity of the spin-phonon interaction and the potential in giving a universal form of the relaxation rate from rotation-only vibration modes, some studies have been done with encouraging results. Specifically, by using the rotational spin-phonon Hamiltonian Dohm and Fulde (1975); Bonsall and Melcher (1976), a universal formula of direct transition rate between tunnel-split states in a strong axially anisotropic non-Kramers system described by ZFS in Eq. (1) has been found Chudnovsky et al. (2005). In the same work, the direct process rate formula for a specific condition $\left[\mathcal{H}_{A},S_{z}\right]=0$ was also provided in spite of its questionable applicability. Within the same approach, Raman process was also considered for a specific non-Kramers system with biaxial spin Hamiltonian Calero et al. (2006). However, as far as we are aware, both direct and Raman process for arbitrary spin Hamiltonians and type of doublet(Kramers and non-Kramers) have not yet been considered despite their importance for the understanding of magnetic relaxation in a broad class of systems.

In this work, we derive universal expressions for the rate of spin-phonon relaxation at low temperature in Kramers and non-Kramers systems in the presence of applied magnetic field based on the “rotational” contribution to the spin-phonon interaction. These include direct process and first- and second-order Raman processes. In order to exclude the effect of hyperfine and dipolar interaction, we restrict the consideration to a strongly diluted systems in which the decoherence is substantially weakened and the system is therefore in the coherent relaxation regime. The article is divided into seven sections. In Section II, spin-phonon relaxation in SMMs at low temperature and the rotational spin-phonon Hamiltonian are introduced; general expressions for the direct and Raman process are also developed. Detailed form of the corresponding relaxation rates for Kramers and non-Kramers system are respectively presented in Section III and IV. Section V is devoted to the comparison between the direct process and Raman processes in different conditions. Section VI shows an example of application of obtained expressions to a real system. A summary of the novel expressions for the relaxation rates and the discussion of their application is given in the last section.

II Spin-phonon relaxation at low temperature and “rotational” spin-phonon Hamiltonian

The spin-phonon relaxation in a SMM is governed by the Redfield equation which in the interaction picture isGaranin (2011),

[TABLE]

where $\rho_{\alpha\beta}$ is the density matrix written in the eigenstates of the ZFS of the ground $S$ term or of CF of the atomic $J$ multiplet and $R_{\alpha\beta,\alpha^{\prime}\beta^{\prime}}$ are matrix elements of the relaxation Redfield super-operator. Due to the complexity of this equation, some approximations are used, the most known being the secular approximation. It is based on the fact that the time change of the density matrix element due to the relaxation is slow compared to the frequency difference $\left|\omega_{\alpha\beta}-\omega_{\alpha^{\prime}\beta^{\prime}}\right|$. As a consequence, the corresponding fast terms on the right-hand side of the equation can be averaged out and ignored. At low temperature, when only the ground doublet eigenstates $\ket{\pm}$ are populated, the off-diagonal density matrix elements $\rho_{+-}$ $\left(\rho_{-+}\right)$ will decay with the rate $\tilde{\Gamma}_{+-}\equiv-R_{+-,+-}=1/2\left(R_{++,--}+R_{--,++}\right)\equiv 1/2\left(\Gamma_{+-}+\Gamma_{-+}\right)$ since the other terms on the right-hand side (RHS) of Eq. (3) oscillate rapidly with the frequency $\left|\omega_{+-}\right|$ or 2$\left|\omega_{+-}\right|$ and, therefore, can be neglected within the secular approximation if $\tilde{\Gamma}_{+-}\ll\left|\omega_{+-}\right|$. For the diagonal density matrix elements $\rho_{++}$ $\left(\rho_{--}\right)$, whereas the left-hand side varies with the rate $\Gamma_{\mathrm{relax}}=\Gamma_{+-}+\Gamma_{-+}$ 222The derivation is given in Sec. IV.E2 of Ref. Garanin, 2011, the terms containing off-diagonal density matrix elements in the RHS oscillate with frequency $\left|\omega_{+-}\right|$ and the condition for ignoring these fast terms turns out to be similar, $\Gamma_{\mathrm{relax}}\ll\left|\omega_{+-}\right|$ , which is thus a joint condition for applying the secular approximation in this regime. This condition seems to be satisfied for most SMMs at sufficiently low temperature.

Defining the “low temperature” in this way, one should have in mind that it also implies that the excited states are sufficiently high that the relaxation proceeds within the ground doublet only. Since all excited doublets are excluded at low temperature, the Orbach relaxation process is ruled out. As evident from the expression $\Gamma_{\mathrm{relax}}=\Gamma_{+-}+\Gamma_{-+}$, the relaxation of the magnetic system now comes only from the direct and Raman process between two eigenstates of the ground doublet, i.e. between its two tunnel-split states. This is opposite to the incoherent quantum tunneling regimeLeuenberger and Loss (2000); Garanin and Chudnovsky (1997); Gatteschi et al. (2006); Abragam and Bleaney (1970) which arises when electron-phonon escape rate to excited states considerably exceeds $\Gamma_{\mathrm{relax}}$. This problem is appropriately described in the natural basis of the two states of the ground doublet, which are $\ket{\pm m}=\frac{1}{\sqrt{2}}\left(\ket{+}\pm\ket{-}\right)|_{\bm{\mathrm{H}}=0}$, and generally requires going beyond the simple secular approximation Ho and Chibotaru (2017). Lowering the temperature, the electron-phonon escape rate is suppressed, however, incoherent quantum tunneling regime can still be maintained due to interaction with nuclear spins and through spin dipolar interaction between magnetic molecules in the lattice Vijayaraghavan and Garg (2009); Vijayaraghavan (2011); Prokof’ev and Stamp (2000); Garanin et al. (2000); Sinitsyn and Prokof’ev (2003), provided the tunneling gap is sufficiently small. Since we do not introduce any restriction for its value, hereinafter we keep out the incoherent quantum tunneling effects caused by coupling with nuclear spins (hyperfine field) or other magnetic molecules (dipolar field), by considering strongly diluted systems. Hence, the calculation of the relaxation rate of SMM systems reduces to calculating the relaxation rate between two tunnel-split states of the ground doublet, which consists of the direct and Raman processes. To this end, a knowledge of spin-phonon Hamiltonian is required.

In the following, the treatment is done explicitly for $S$-complexes. However, given their generality, the results are straightforwardly applied also to systems characterized by total angular momentum $J$ as well. That is to say, the results are valid for both transition metal and lanthanide-based magnetic complexes.

As mentioned, a magnetic molecule interacts with the lattice through molecular degrees of freedom including both distortions and rotations. For a rigid magnetic molecule or cluster, the “rotational” contribution, which changes the direction of the main anisotropy axis, is apparently the most efficient. In the presence of an external magnetic field, the total Hamiltonian for the system interacting in this way with the phonon bath can be written as

[TABLE]

where $R=e^{-i\bm{\mathrm{S}}\cdot\bm{\mathrm{\delta\varphi}}}$ is the rotation operator corresponding to a small rotation angle $\bm{\mathrm{\delta\varphi}}$ of the entire complex Chudnovsky et al. (2005); Garanin (2011), $\mathcal{H}_{A}$ is the anisotropic crystal (ligand) field Eq. (1) or (2), $\mathcal{H}_{\mathrm{ph}}$ is the phonon bath Hamiltonian, and $\mathcal{H}_{Z}=g\mu_{B}\bm{\mathrm{H}}\cdot\bm{\mathrm{S}}$ ($\mathcal{H}_{Z}=g_{J}\mu_{B}\bm{\mathrm{H}}\cdot\bm{\mathrm{J}}$ in the case of lanthanides) is the Zeeman Hamiltonian corresponding to the external magnetic field $\bm{\mathrm{H}}$. Expanding Eq. (4) up to the second order in $\bm{\delta\varphi}$ yields the spin-phonon Hamiltonian Garanin (2011),

[TABLE]

As discussed in the introduction, we further consider the contribution of acoustic phonons only as being the most efficient in direct and Raman process at relatively low temperaturesChudnovsky et al. (2005). The angle $\bm{\delta\varphi}$ then can be quantized using the canonical phonon quantization of the lattice displacement vector $\bm{\mathrm{u}}\left(\bm{\mathrm{r}}\right)$ Leuenberger and Loss (2000); Garanin (2011),

[TABLE]

Here $\rho$ is the mass density, $V$ is the volume of the crystal, $\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ is the phonon frequency, and $\bm{\mathrm{e}}_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ is the polarization unit vector corresponding to one longitudinal and two transverse polarizations, $\lambda=t_{1},t_{2},l$. As can be seen from the above expression for $\delta\bm{\mathrm{\varphi}}$, the longitudinal mode does not contribute to the spin-phonon Hamiltonian due to the fact that $\bm{\mathrm{k}}\parallel\bm{\mathrm{e}}_{\bm{\mathrm{k}}l}$. Since we consider a low temperature, the Debye model’s dispersion relation $\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}=v_{\lambda}k$ will be used throughout this work.

Denoting the states of the thermal bath by the phonon population numbers $n_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ for each phonon mode $\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}$, the direct and Raman transition rates between the tunnel-split eigenstate $\ket{+}$ and$\ket{-}$ of the ground doublet can be calculated from the Fermi’s golden rule Cohen-Tannoudji et al. (1977),

[TABLE]

where $\Omega\equiv\varepsilon_{+}-\varepsilon_{-}$ and

[TABLE]

Here $\ket{\xi}$ labels eigenstates of the spin $S$. The transition matrix element $M^{\left(1\right)}$ characterizes the direct process where owing to the interaction Hamiltonian $\mathcal{H}^{\left(1\right)}$, the spin system emits one phonon to transit from $\ket{+}$ to $\ket{-}$, while $M^{\left(2\right)}$ characterizes the first-order Raman process where the spin system under the effect of $\mathcal{H}^{\left(2\right)}$ virtually emits a phonon $\hbar\omega_{\bm{\mathrm{q}}}$ and concomitantly absorbs another one with energy $\hbar\omega_{\bm{\mathrm{k}}}$, their energy difference spanning the states $\ket{+}$ and $\ket{-}$ (Fig. 1); $M^{\left(1+1\right)}$ characterizes the second-order Raman process with the same physical interpretation as the first-order Raman process but deduced from the second order of perturbation of the spin-phonon Hamiltonian $\mathcal{H}^{\left(1\right)}$ (Fig. 1).

In order to determine the relaxation rate in Eq. (9) and (10)), one should calculate the matrix elements $M^{\left(i\right)}$.

Direct process:

From Eq. (6) and (11), we easily reach

[TABLE]

where

[TABLE]

Substituting $M^{\left(1\right)}$ into Eq. (9) of $\Gamma_{-+}^{\mathrm{dr}}$ and averaging over the number of phonon numbers $n_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ results in

[TABLE]

Using the property $\left[\bm{\mathrm{k}}\times\bm{\mathrm{e}}_{\bm{\mathrm{k}}t}\right]=\pm k\bm{\mathrm{e}}_{\bm{\mathrm{k}}t^{\prime}}$ where $t=t_{1},t_{2}$ and $t^{\prime}=t_{2},t_{1}$ respectively and the Debye dispersion relation for the acoustic phonons, then replacing $\sum_{\bm{\mathrm{k}}}$ by $\frac{V}{8\pi^{3}}\int d^{3}k$, one comes up to a simple expression for $\Gamma_{-+}^{\mathrm{dr}}$ (see Appendix A for detailed calculation):

[TABLE]

Here $v_{t_{1}}$, $v_{t_{2}}$ are respectively the phonon velocity for two transverse polarizations, $t_{1}$ and $t_{2}$. The expression for $\bm{\mathrm{\Xi}}=\braket{-}{\mathcal{H}^{\left(1\right)}}{+}$ is different for Kramers and non-Kramers doublet.

Raman process:

From Eq. (12) for $M^{\left(2\right)}$, we obtain

[TABLE]

where

[TABLE]

Similar derivation for $M^{\left(1+1\right)}$ yields

[TABLE]

During the derivation, the intermediate states $\ket{n_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}-1,n_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}}$ and $\ket{n_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}},n_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}+1}$ of the bath, which correspond to the first and second term in the above equation, have been substituted into Eq. (13).

Since $M^{\left(1+1\right)}$ connects the ground doublet eigenstate $\ket{+}$ to $\ket{-}$ via the states $\ket{\xi}$, it is not necessary to take into account the admixture of other excited doublets caused by the Zeeman interaction. Denoting this non-admixing part of the eigenstates $\ket{\pm}$ by $\ket{p_{1}\left(m_{1}\right)}$ (short notation for “plus” and “minus”) and the non-admixing part of $\ket{\xi_{m}}$ corresponding to $m^{\mathrm{th}}$ doublet by $\ket{x_{m}}$, $x=p$ or $m$, we have,

[TABLE]

where $\delta M_{\mathrm{sp}}^{\alpha\beta}$ results from the excited doublets admixing part of the eigenstates $\ket{\pm}$ in the expression of $M^{\left(2\right)}$. Its detailed form will be shown separately for Kramers and non-Kramers system.

Finally, we set two approximations further used for calculation of the Raman relaxation rate. First, since we are working in low temperature regime, it is assumed that $\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}\ll\hbar\omega_{\xi 1}$, $\ket{\xi}\neq\ket{\pm}$. Second, due to the nature of the second-order Raman process, most of the phonon frequencies involved in the expression for$M^{\left(1+1\right)}$ are supposed to be much larger than the energy gap between $\ket{+}$ and $\ket{-}$, i.e. $\Omega\ll\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$. We thus have the combined assumption,

[TABLE]

III Spin-phonon relaxation in Kramers doublets

As mentioned above, at low temperature only the ground doublet is populated and the relaxation rate of the system reduces to the one between the two tunnel-split states. For Kramers system, without an applied magnetic field, the direct and first-order Raman process disappear in virtue of the Kramers’ theorem Abragam and Bleaney (1970). The latter is violated in the presence of an external magnetic field through the admixture of excited Kramers doublets via Zeeman interaction Abragam and Bleaney (1970). This process thus requires an enlargement of the working Hilbert space to include all excited doublets of $S$ even though these are not thermally populated. In this space, the electronic Hamiltonian in zero external magnetic field can be written as

[TABLE]

where $\ket{\bar{m}^{\left(0\right)}}=\theta\ket{m^{\left(0\right)}}$ is the Kramers conjugate states of the eigenstate $\ket{m^{\left(0\right)}}$, $\theta$ is the time-reversal operator, and the index $m=1\ldots\left(S+1/2\right)$ indicates the doublets. Since we are interested in the ground doublet, a choice of the reference frame’s axes along three main magnetic axes of this doublet will apparently be the most appropriate. Further, the main magnetic $z$-axis of the pseudospinChibotaru and Ungur (2012) $\tilde{S}=1/2$ corresponding to the ground doublet will be chosen as the quantization axis. The submatrix of the electronic Hamiltonian in a magnetic field $\bm{\mathrm{H}}=(H_{x},H_{y},H_{z})$ then is

[TABLE]

where

[TABLE]

Here $g_{x},g_{y},g_{z}$ are the principal values of the $g$-tensor of the ground doublet. For convenience, we set $\mu_{B}=1$ and $\varepsilon_{1}=0$.

To make this Hamiltonian submatrix real, the basis vectors of the ground doublet are rotated by an angle $\varphi/2$ while other basis vectors are kept unchanged:

[TABLE]

Hereinafter, we will use $M$ to denote all basis vectors and $m$ to denote states corresponding to $m^{\mathrm{th}}$ doublet.

In this new basis, the submatrix $\left(\mathcal{H}_{A}+\mathcal{H}_{Z}\right)_{\left\{\ket{1},\ket{\bar{1}}\right\}}$ becomes real,

[TABLE]

In order to find the spin-phonon coupling matrix element of the Kramers doublet in the presence of an external magnetic field, we need to pass to the new eigenstates of the spin system. These are given by the diagonalization of the submatrix $\left(\mathcal{H}_{A}+\mathcal{H}_{Z}\right)_{\left\{\ket{1},\ket{\bar{1}}\right\}}$:

[TABLE]

where $\vartheta$ is defined as

[TABLE]

As can be seen, $\ket{p_{1}}$ and $\ket{m_{1}}$ are still Kramers conjugates states. Accordingly, without admixture of excited doublets, the matrix element of the spin-phonon Hamiltonian, which is a time-even operator, between them is still zero. Applying the perturbation theory w.r.t Zeeman Hamiltonian $\mathcal{H}_{Z}$, we get new eigenstates (see Appendix B for details):

[TABLE]

For the later use, we also give here the expressions for $\ket{p_{m}}$ and $\ket{m_{m}}$:

[TABLE]

For the calculation of the direct and Raman transition rates, we need to find the matrices of spin operators $\bm{\mathrm{S}}$ (or $\bm{\mathrm{J}}$ for lanthanide complexes). Notice that in the original sub-basis $\left\{\ket{1^{\left(0\right)}},\ket{\bar{1}^{\left(0\right)}}\right\}$, the matrix of the magnetic moment $\mu_{\alpha}$ can be calculated by either $\mu_{\alpha}^{\left(0\right)}=-gS_{\alpha}^{\left(0\right)}$ or $\mu_{\alpha}^{\left(0\right)}=-g_{\alpha}\tilde{S}_{\alpha}^{\left(0\right)}=-\frac{1}{2}g_{\alpha}\sigma_{\alpha},$ where $g$ is the isotropic electronic $g$-factor of the whole spin multiplet $S$ and $\sigma_{\alpha}$ $(\alpha=x,y,z)$ are Pauli matrices (in the case of $J$-complexes one uses for the first equation $\mu_{\alpha}^{\left(0\right)}=-g_{J}J_{\alpha}^{\left(0\right)},$ where $g_{J}$ is the Landé $g$-factor), we come to a simple formula relating the real-spin $S_{\alpha}$ and pseudospin $\tilde{S}_{\alpha}$ in the original sub-basis $\left\{\ket{1^{\left(0\right)}},\ket{\bar{1}^{\left(0\right)}}\right\}$:

[TABLE]

A change from sub-basis $\left\{\ket{1^{\left(0\right)}},\ket{\bar{1}^{\left(0\right)}}\right\}$ to $\left\{\ket{1},\ket{\bar{1}}\right\}$, therefore, results in a rewriting of the $\bm{\mathrm{S}}$ matrices

[TABLE]

Equipped with these preliminary results, we now proceed to the calculation of the direct and Raman transition rates.

III.1 The direct process

For the direct transition rate $\Gamma_{-+}^{\mathrm{dr}}$, the value of $\bm{\mathrm{\Xi}}=\braket{-}{\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}}{+}$ is required. Substituting the eigenstates $\ket{\pm}$ from Eq. (36,37) into $\bm{\mathrm{\Xi}}$, then expanding it to the 1st order of perturbation theory, yields

[TABLE]

where the property $\braket{m_{1}}{\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}}{p_{1}}=0$ ($\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}$ is a time-even operator) has been used (see Appendix C).

Substituting the expressions for $\ket{p_{1}}$ and $\ket{m_{1}}$ from Eq. (33) and (34) into the above equation and making use of the time reversal symmetry in the matrix elements of time-even/time-odd operators (see Appendix C), the equation is simplified to

[TABLE]

The explicit form of $\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}=i\left[\mathcal{H}_{A},\bm{\mathrm{S}}\right]$ gives a key relation,

[TABLE]

Inserting this into Eq. (45) results in

[TABLE]

The expression for the vector $\bm{\mathrm{z}}_{1}$ can be simplified by using the completeness relation, properties of the time-even operator Abragam and Bleaney (1970), and the relation $\left[\mathcal{H}_{Z},\bm{\mathrm{S}}\right]=ig\left(\bm{\mathrm{S}}\times\bm{\mathrm{H}}\right)$:

[TABLE]

From the explicit form of $\bm{\mathrm{S}}$ matrices in Sec. II, we have $\bm{\mathrm{S}}_{11}=-\bm{\mathrm{S}}_{\bar{1}\bar{1}}=\left(0,0,g_{z}\right)/2g$ and $\bm{\mathrm{S}}_{1\bar{1}}=\bm{\mathrm{S}}_{\bar{1}1}^{*}=\left(g_{x}e^{i\varphi},-ig_{y}e^{i\varphi},0\right)/2g$. This results in

[TABLE]

where

[TABLE]

Similarly, we have for $\bm{\mathrm{z}}_{2}$ and $\bm{\mathrm{z}}_{3}$,

[TABLE]

Inserting these vectors into Eq. (47) and then into Eq. (18), we obtain the final expression for the direct transition rate between eigenstates of a Kramers doublet,

[TABLE]

A remarkable feature of this expression is that it does not involve explicitly any parameters characterizing the excited doublets despite the fact that their admixture is indispensable for the existence of the transition. This is due to the special form of the rotational spin-phonon Hamiltonian mentioned earlier and the implicit incorporation of the contributions of the admixed excited doublets in the ground doublet’s $g$-tensor, which thus allows a very compact and handy form of $\Gamma_{-+}^{\mathrm{dr}}$.

Eq. (55) shows that the direct transition rate between the eigenstates of the doublet involves three contributions. The first contribution depends solely on the energy bias $W$ between two eigenvectors of the pseudospin operator $\tilde{S}_{z}$ and hence may be considered as the direct process between two Kramers conjugate states $\ket{1}$ and $\ket{\bar{1}}$ of the ground doublets. . The third contribution only depends on the tunneling splitting gap $\Delta$ and can be regarded as the direct process between two tunneling states $\ket{\pm}$ in the absence of the bias. Meanwhile, the second contribution appears as a term mixing these two.

The $\Omega$-dependence of the relaxation in (55) looks surprisingly different from the conventional expression for a direct transition rate between two Kramers components Abragam and Bleaney (1970). However, in the case of $\Delta\ll W$, only the first contribution survives, we have $\Omega\approx W$ and Eq. (55) acquires a conventional prefactor $\Omega^{5}$. Thus the derived direct transition rate, besides being of a universal form (i.e., not containing explicitly the parameters of spin-phonon coupling Chudnovsky et al. (2005)), is also a generalization of a previous expression Abragam and Bleaney (1970) in that it contains distinct contributions of the bias ($W$) and tunneling gap ($\Delta$), which thus cannot be described as depending on only one direct gap $\Omega=\sqrt{\Delta^{2}+W^{2}}$ as was thought before Abragam and Bleaney (1970).

Since $W,\Delta\propto H$, it is obvious that the derived direct transition rate is proportional to $H^{4}T$ when $\Omega\ll k_{\mathrm{B}}T$. This is thus in agreement with the rough approximation for the field and temperature dependencies of the direct process in Kramers doublet in Ref. [Abragam and Bleaney, 1970]. However, Eqs. (26) and (27) show that all contributions to $\Gamma_{-+}^{\mathrm{dr}}$ in (55) depend also on the direction of the applied magnetic field.

In the case $g_{x}=g_{y}\equiv g_{\perp}$ and $\Omega\ll k_{\mathrm{B}}T$, the transition rate is explicitly independent of the magnetic field orientation and proportional to the square of $g_{\perp}$,

[TABLE]

This physically means that the more axially anisotropic the system is, the smaller is the direct process relaxation rate.

On the other hand, for the near-isotropic system such that $g_{x}\approx g_{y}\approx g_{z}\approx g$, resulting in $g_{i}\ll g_{\alpha}$ ($i=1,2,3$ and $\alpha=x,y,z$), $\Gamma_{-+}^{\mathrm{dr}}$ being proportional to $g_{i}$ asymptotically approaches zero in the limit of totally isotropic spin system.

III.2 The Raman process

For the evaluation of $\Gamma_{-+}^{\mathrm{Raman}}$, we need to calculate $M^{\left(2\right)}+M^{\left(1+1\right)}$. Since $\ket{m_{1}}=\ket{\bar{p}_{1}}$, Eq. (34), the first term of $M_{\mathrm{sp}}^{\alpha\beta}$ in Eq. (22) for $M^{\left(2\right)}+M^{\left(1+1\right)}$ is zero due to time-reversal symmetry. Meanwhile, the third term can be decomposed into two components.The first involves states of the ground doublet only and therefore is zero due to time-reversal symmetry as well. The second component of the third term, dubbed $M_{\mathrm{sp2}}^{\alpha\beta}$, can be simplified using the assumption (23):

[TABLE]

where, like in the case of direct process, Eq. (44), $\hbar\omega_{1m}$ denotes the excitation energy between the ground and $m^{\mathrm{th}}$ excited Kramers doublets in the absence of applied field. Here, we made the approximation $\delta\left(\Omega+\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}-\hbar\omega_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}\right)\approx\delta\left(\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}-\hbar\omega_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}\right)$ and thus replaced $\omega_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}$ by $\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$.

It is easy to show that the first term of $M_{\mathrm{sp2}}^{\alpha\beta}$ is zero. Indeed, for any $m^{\mathrm{th}}$ doublet (in which $\ket{m_{m}}=\ket{\bar{p}_{m}}$), the numerator of the first term becomes

[TABLE]

The second term on the r.h.s of Eq. (57) can be easily found by using the property $\mathcal{H}_{A}\ket{x_{m}}=\varepsilon_{m}\ket{x_{m}}$ and $\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}=i\left[\mathcal{H}_{A},\bm{\mathrm{S}}\right]$,

[TABLE]

where $\varepsilon_{\alpha\beta\gamma}$ is the Levi-Civita tensor.

The second term $\delta M_{\mathrm{sp}}^{\alpha\beta}$, in the expression for $M^{\left(2\right)}+M^{\left(1+1\right)}$, Eq. (22), has the explicit form:

[TABLE]

where again the denominator corresponds to the absence of applied field.

It can be easily seen that $\delta M_{\mathrm{sp}}^{\alpha\beta}$ is of the order $\mathcal{O}\left(\mathcal{H}^{Z}S^{2}\right)$, whereas, from Eq. (59), $M_{\mathrm{sp2}}^{\alpha\beta}$ is of the order $\mathcal{O}\left(\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}S^{2}\right)$. Since the Zeeman interaction is considered to be much smaller than $\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ according to the condition (23), $\delta M_{\mathrm{sp}}^{\alpha\beta}$ can be safely neglected in comparison to $M_{\mathrm{sp2}}^{\alpha\beta}$. As a consequence,

[TABLE]

Averaging over the phonon number $n_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$, $n_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}$ and summing up over all phonon modes (see Appendix A for details) as in the calculation of the direct process transition rate, we finally obtain,

[TABLE]

where $\omega_{D}$ is the Debye frequency and

[TABLE]

The value of $Q$ can be calculated straightaway from the explicit form of $\ket{p_{1}}$ Eq. (33), $\ket{m_{1}}$ Eq. (34), and $S_{\alpha}$ Eq. (41-43),

[TABLE]

where $\sin\varphi=\Delta_{y}/\Delta$, $\cos\varphi=\Delta_{x}/\Delta$, $\sin\vartheta=\Delta/\Omega$, and $\cos\vartheta=W/\Omega$.

This yields,

[TABLE]

As can be seen, the Raman process transition rate varies as $T^{9}$ as expected. Surprisingly, the Raman process transition rate only shows a dependence on the orientation of the applied field but not on its magnitude. An interesting consequence of this property is that in experiments with powder sample of a rigid Kramers molecules, the Raman process transition rate is independent of the applied magnetic field.

Interestingly, for the near-isotropic systems with $g_{x}\approx g_{y}\approx g_{z}\approx g$, the Raman process transition rate derived above yields a small value due to its proportionality to $g_{i}^{2}$, which are expected to be very small according to Eqs. (50-52). This is in line with the situation for the direct process.

If the Debye frequency $\omega_{D}$ is much larger than $T$, $I_{8}\approx 8!$ and

[TABLE]

In the case the energy bias $W$ exceeds the tunneling splitting $\Delta$, a simpler expression can be obtained:

[TABLE]

On the other hand, if $g_{x}=g_{y}\equiv g_{\perp}$ and accordingly $g_{1}=g_{2}\equiv g_{12}$, Eq. (65) becomes

[TABLE]

As can be seen, at variance to a similar situation in the case of the direct process, Eq. (56), $\Gamma_{-+}^{\mathrm{Raman}}$ still depends on the direction of applied field.

IV Spin-phonon relaxation in non-Kramers doublets

A similar treatment can also be applied to non-Kramers systems with a little modification. Particularly, for non-Kramers (quasi) doublets, we introduce a new basis $\ket{m^{*}}$ and $\ket{m^{\prime*}}$ of $m^{\mathrm{th}}$ doublet:

[TABLE]

where $\ket{\pm_{m}^{\left(0\right)}}$ are eigenstates of the anisotropic Hamiltonian $\mathcal{H}_{A}$ corresponding to the $m^{\mathrm{th}}$ doublet, and related to the corresponding real-valued eigenstates, denoted $\ket{\pm_{m,r}^{\left(0\right)}}$, as

[TABLE]

Obviously,

[TABLE]

In this basis, the anisotropic Hamiltonian $\mathcal{H}_{A}$ is of the form:

[TABLE]

where $\Delta_{m}$ is the tunneling gap of the $m^{\mathrm{th}}$ doublet. Hereinafter, we choose $\varepsilon_{1}=0$ and denote all quantities corresponding to the ground doublet ($m=1$) without the index 1, e.g. $\Delta\equiv\Delta_{1}$ or $\ket{\pm_{1}}\equiv\ket{\pm}$.

Choosing the $z$-axis to lie along the main magnetic axis of the ground doublet, the spin Hamiltonian in the sub-basis $\left\{\ket{1^{*}},\ket{\bar{1}{}^{*}}\right\}$ of the doublet under an applied magnetic field becomes

[TABLE]

where $g_{z}$ is the principal value along the main magnetic $z$-axis of the ground doublet’s $g$-tensor Griffith (1963). For non-Kramers doublet, $g_{x}=g_{y}=0$ Griffith (1963).

New eigenstates of the ground doublet can be found by first diagonalizing the submatrix $\left(\mathcal{H}_{A}+\mathcal{H}_{Z}\right)_{\left\{\ket{1^{*}},\ket{\bar{1}^{*}}\right\}}$ then using the perturbation theory,

[TABLE]

Here $\hbar\omega_{MN}\equiv\varepsilon_{M}-\varepsilon_{N}$. Further, we will make use of the expressions for $\ket{p_{m}}$ and $\ket{m_{m}}$:

[TABLE]

Similar to the Kramers system, we also need to know the matrix representation of the spin operators $S_{\alpha}$ in the sub-basis $\left\{\ket{1^{*}},\ket{\bar{1}^{*}}\right\}$. The former choices of the $z$-axis and the basis yield Griffith (1963)

[TABLE]

IV.1 The direct process

As is clear from Eq. (18), the only difference in the expression for $\Gamma_{-+}^{\mathrm{dr}}$ for non-Kramers system compared to Kramers one is the expression for $\bm{\mathrm{\Xi}}$. Substituting $\ket{\pm}$ into $\bm{\mathrm{\Xi}}=\braket{-}{\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}}{+}$, we have

[TABLE]

Further derivation is practically identical to the Kramers’ case. The key difference here is that the term $\braket{m_{1}}{\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}}{p_{1}}$ is non-zero for non-Kramers doublet. Taking into account that for non-Kramers system $\theta^{2}=1$ and the matrix elements from time-even and time-odd operators possess the properties presented in Appendix C, we come to the following result:

[TABLE]

$\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}=i\left[\mathcal{H}_{A},\bm{\mathrm{S}}\right]$ can be expanded as follows:

[TABLE]

Since $\Delta_{M}$ for all $M$ is much smaller than the energy gap between any two doublets, we have

[TABLE]

Using the completeness relation and $\left[\mathcal{H}_{Z},\bm{\mathrm{S}}\right]=ig\left(\bm{\mathrm{S}}\times\bm{\mathrm{H}}\right)$, we obtain

[TABLE]

The term $\bm{\mathrm{\mathcal{H}}}_{m_{1}p_{1}}^{\left(1\right)}$ can be found by utilizing the easy-to-prove property,

[TABLE]

then simplifying,

[TABLE]

Hence,

[TABLE]

and

[TABLE]

It is easy to see that when $\Omega\ll k_{\mathrm{B}}T$ and the energy bias $W$ dominates over the tunneling splitting $\Delta$, the direct transition rate $\Gamma_{-+}^{\mathrm{dr}}$ is proportional to $H^{2}T$ as expected for non-Kramers systems Abragam and Bleaney (1970).

Interestingly, since the rate is proportional to the square of the tunneling splitting $\Delta$, a reduction in the latter will substantially decrease the direct process transition rate. That is to say, the more anisotropic the system, the slower the direct process transition rate. This is not directly evident in the standard expression for the direct relaxation rate Abragam and Bleaney (1970).

In a highly anisotropic system where $W=g_{z}H_{z}\approx gH_{z}\left(m-m^{\prime}\right)$ where $m$ ($m^{\prime}$) is the magnetic quantum numbers corresponding to the ground doublet, the above expression becomes

[TABLE]

which is similar to the result in Ref. [Chudnovsky et al., 2005] if $v_{t_{1}}=v_{t_{2}}$. Compared to it, Eq. 93 is more general because it does not require extreme uniaxial anisotropy and is applicable not only to $S$ systems (transition metal complexes) but also to $J$ systems (lanthanide complexes) as well.

IV.2 The Raman process

For non-Kramers system, the term $M_{\mathrm{sp}}^{\alpha\beta}$ in Eq. (22) can be expanded as,

[TABLE]

Using the property (90) and $\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}=i\left[\mathcal{H}_{A},\bm{\mathrm{S}}\right]$, it can be easily proved that $\mathcal{H}_{p_{1}p_{1}}^{\left(1\right)\alpha}=\mathcal{H}_{m_{1}m_{1}}^{\left(1\right)\alpha}=0$ for all $\alpha$. Consequently,

[TABLE]

Similar to the calculation of $M_{\mathrm{sp1}}^{\alpha\beta}$, noticing that all $\Delta_{i}$ can be neglected as being much smaller than $\hbar\omega_{1m}$, $M_{\mathrm{sp2}}^{\alpha\beta}$ becomes 333Even if we keep all terms containing $\Delta_{i}$, the order of magnitude of these terms, which is $\mathcal{O}\left(\Delta_{i}S^{2}\right)$, is still smaller than one of $M_{\mathrm{sp4}}^{\alpha\beta}$, which is $\mathcal{O}\left(\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}S^{2}\right)$ as we see later. Hence, by any means, these terms can be safely neglected.

[TABLE]

Likewise, $M_{\mathrm{sp3}}^{\alpha\beta}$ can be simplified into:

[TABLE]

Here, we have ignored all terms containing $\Delta_{i}$ as before.

$M_{\mathrm{sp4}}^{\alpha\beta}$ can be determined without difficulty,

[TABLE]

The term $\delta M_{\mathrm{sp}}^{\alpha\beta}$ for non-Kramers system has the explicit form:

[TABLE]

As in the Kramers’ case, this term is of the magnitude order $\mathcal{O}\left(\mathcal{H}^{Z}S^{2}\right)$, which is smaller than the order of magnitude of $M_{\mathrm{sp4}}^{\alpha\beta}$, $\mathcal{O}\left(\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}S^{2}\right)$, since characteristic $\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ is much larger than the Zeeman interaction. As a consequence, $\delta M_{\mathrm{sp}}^{\alpha\beta}$ can be ignored. Then $\Gamma_{-+}^{\mathrm{Raman}}$ for non-Kramers system acquires a similar form, (61), with the Kramers case. Performing similar calculations as for the Kramers’ case results in the same Eq. (62) but with a different value of $Q\equiv\sum_{\alpha,\beta}\left|Q_{\alpha\beta}\right|^{2}$,

[TABLE]

resulting in

[TABLE]

If the Debye frequency $\omega_{D}$ is much larger than $T$, we have

[TABLE]

In a highly anisotropic system where $W=g_{z}H_{z}\approx gH_{z}\left(m-m^{\prime}\right)$, the above expression becomes

[TABLE]

The results obtained here differ the previous treatment of the Raman process in non-Kramers doublets within the same rotational approximation for the spin-phonon coupling in several important aspects. Contrary to the present treatment, Ref. [Calero et al., 2006] is limited to a specific axially-anisotropic spin Hamiltonian $\mathcal{H}_{A}$, also called biaxial spin Hamiltonian. It is also important to note that the present theory predicts $T^{9}$ dependence of the Raman relaxation rate for non-Kramers doublets (within rotational spin-phonon approximation), while the previous treatment shows a $T^{11}$ temperature dependence Calero et al. (2006). The discrepancy comes from the zero value of the last square bracket in the expression (95) for $M_{\mathrm{sp}}^{\alpha\beta}$, found in [Calero et al., 2006], but not supported by the present calculations.

V Direct process versus Raman process in magnetic molecules

The derived expressions for direct process and Raman process transition rates permit us to assess the relative importance of these processes in the typical case of $\Omega\ll k_{\mathrm{B}}T$ and $\hbar\omega_{D}\gg k_{\mathrm{B}}T$. This analysis is valid for cases when rotational spin-phonon interaction can be considered dominant

V.1 Kramers SMMs

Since there are many parameters in the expressions for $\Gamma_{-+}^{\mathrm{dr}}$ Eq. (55) and $\Gamma_{-+}^{\mathrm{Raman}}$ Eq. (65), we will take as example a typical lanthanide SMM with $g_{J}=4/3$ $(\mathrm{Dy}^{3+})$, $g_{x}=g_{y}=g_{\perp}=0.01$, $g_{z}\approx 20$ under the supposition of equal $H_{\alpha}$ ($\alpha=x,y,z$). Taking for the common factor in Eqs. (55) and (65), $E_{t}\equiv\left[2\hbar^{3}\rho v_{t_{1}}^{5}v_{t_{2}}^{5}/\left(v_{t_{1}}^{5}+v_{t_{2}}^{5}\right)\right]^{1/4}$, a typical value $E_{t}\approx 100\,\mathrm{K}$, $H$ is measured in kOe, we obtain

[TABLE]

Under an applied field of 1 kOe, the Raman transition rate becomes smaller than the direct one only at $T<2.33$ K. Particularly, at $T=5$ K, the Raman transition rate is $1.28\,\mathrm{s^{-1}}$ whereas the direct one is three orders of magnitudes smaller, $0.003\,\mathrm{s^{-1}}$. Fig. 2a shows the dependence of the relaxation times $\tau_{\mathrm{direct}}=\left(2\Gamma_{-+}^{\mathrm{dr}}\right)^{-1}$, $\tau_{\mathrm{Raman}}=\left(2\Gamma_{-+}^{\mathrm{Raman}}\right)^{-1}$, and $\tau_{\mathrm{total}}=\left[2\left(\Gamma_{-+}^{\mathrm{dr}}+\Gamma_{-+}^{\mathrm{Raman}}\right)\right]^{-1}$ as a function of $T^{-1}$ for an applied field $H=1$ kOe, Fig. 2b presents the field dependence of the ratio $\Gamma_{-+}^{\mathrm{Raman}}/\Gamma_{-+}^{\mathrm{dr}}$ for several values of $T$, and Fig. 2c shows the field dependence of these relaxation times at $T=5$ K.

Several features can be inferred from Fig. 2. First, it can be seen (Fig. 2a) that the switch from one process to another occurs sharply in a small domain of a few Kelvin. This is basically due to the fact that the direct process depends linearly on $T$ while the Raman one shows a $T^{9}$ dependence. Additionally, as expected, the direct process starts gaining control over the Raman process when the field increases (Fig. 2b and 2c). However, for a Kramers system the gain is rather slow. Moreover, the transition from dominant Raman to dominant direct process occurs without passing through an extremum of the total relaxation time/rate, a situation quite opposite to the case of non-Kramers SMMs (see below).

V.2 Non-Kramers SMMs

For non-Kramers systems, the number of parameters needed to describe the relaxation is substantially reduced. In fact, in most of cases, we can approximate $\left|\bm{\mathrm{S}}_{1^{*}1^{*}}\right|\approx S$ (or for a lanthanide SMM, $\left|\bm{\mathrm{J}}_{1^{*}1^{*}}\right|\approx J$). Furthermore, since the transverse magnetic field hardly affects the tunneling splitting $\Delta$, only the field projection $H_{z}$ directed along the easy axis of the SMM will affect the relaxation. Then we have

[TABLE]

where $E_{t}$ is defined in the previous subsection.

Accordingly, a $\mathrm{Tb^{3+}}$ SMM with $g_{J}=3/2$, $J=6$, $E_{t}=100$ K, and $\Delta=10^{-2}$ K, will be characterized by the following relaxation rates:

[TABLE]

where $H_{z}$ is in kOe. In zero field, $\Gamma_{-+}^{\mathrm{Raman}}/\Gamma_{-+}^{\mathrm{dr}}\approx\left(T/0.48\right)^{8}$. As can be seen, even with relatively large tunneling splitting gap of $10^{-2}\,\mathrm{K}$, the Raman process still dominates over the direct process in zero field. At $T$=5 K, the direct relaxation rate is still very small, $6.25\times 10^{-5}\mathrm{\,s^{-1}}$, while the Raman process relaxation rate is eight orders of magnitude larger, $\Gamma_{-+}^{\mathrm{Raman}}=8.32\times 10^{3}\mathrm{\,s^{-1}}$. However, in a field of 1 kOe, the situation changes significantly,

[TABLE]

We can see that the direct process dominates at $T<5.31\,\mathrm{K}$. Thus, at $T$=3 K, $\Gamma_{-+}^{\mathrm{dr}}=5.49\times 10^{-1}\,\mathrm{s^{-1}}$ is two orders of magnitude larger than $\Gamma_{-+}^{\mathrm{Raman}}=5.73\times 10^{-3}\,\mathrm{s^{-1}}$. This property actually comes from the radically different behavior of the two processes in applied field where the direct process increases fast with the magnetic field but the Raman process is significantly reduced out of resonance.

In analogy to Fig. 2, Fig. 3 shows the field and temperature dependence of the relaxation times $\tau_{\mathrm{direct}}=\left(2\Gamma_{-+}^{\mathrm{dr}}\right)^{-1}$, $\tau_{\mathrm{Raman}}=\left(2\Gamma_{-+}^{\mathrm{Raman}}\right)^{-1}$, and $\tau_{\mathrm{total}}=\left[2\left(\Gamma_{-+}^{\mathrm{dr}}+\Gamma_{-+}^{\mathrm{Raman}}\right)\right]^{-1}$. As seen in Fig. 3a, the switch from one dominant process to another also takes place sharply in a narrow range of temperatures as in the case of Kramers doublets. This also results from the different $T$-dependence of the two relaxation processes. However, in contrast to Kramers systems, the total relaxation time shows an extremum in its $H$-dependence, which is due to the fact that the Raman process transition rate decreases fast with the field. This is not surprising given that the Raman transition rate depends on the ratio $\Delta^{2}/\left(\Delta^{2}+W^{2}\right)$.

Given the typical set of parameters used here, it is likely that for Kramers SMMs, the Raman process will dominate the direct one at temperature $T>5\,\mathrm{K}$ under not too strong magnetic field, $H<2\,\mathrm{kOe}$. For non-Kramers system, it is probable that at $H>1$ kOe and low temperature ($T<5$ K), the direct process takes over the Raman process. However, at temperature $T>10\,\mathrm{K}$ under $H<2\,\mathrm{kOe}$, the Raman process likely suppresses the direct one.

VI Application: relaxation in a $\mathrm{Co^{2+}}$ complex

As an example of applications to real compounds, in Fig. 4a we show the calculated total relaxation rate with the present expressions and the comparison to experimental data for a cobalt(II) complexGómez-Coca et al. (2014) showing slow relaxation 444In cases when optical phonons contribute significantly to the spin-lattice relaxation, the Raman relaxation rate for Kramers system may deviate from the $T^{9}$ dependence Shrivastava1983. Therefore, in order to have an accurate assessment on the applicability of the rotational spin-phonon approximation, only those experiments which satisfies the following criteria are chosen: (1) Raman process shows a $T^{9}$ dependence, which thus excludes the role of optical phonon; (2) a small number of parameters in the expression of the total relaxation rate, which thus prevents over-parameterization; (3) data and fittings are of high-quality. Ref. [Gómez-Coca et al., 2014] satisfies all these criteria.. The cobalt(II) ion possesses a ground state spin $S=3/2$ with an almost isotropic $g$-factor $g\approx 2.5$; the principal values of the ground doublet’s $g$-tensor $g_{x}=2.65$, $g_{y}=6.95$, $g_{z}=1.83$, and the mass density $\rho=1528\,\mathrm{kg/m^{3}}$. Since the experiment was conducted on polycrystalline samples, for a simple account of this in the calculation, the field was oriented at an angle $\theta_{\bm{\mathrm{H}}}=\phi_{\bm{\mathrm{H}}}=\pi/4$ w.r.t anisotropic axes $z$ and $x$. A fitting with Eq. 65 for $\Gamma_{-+}^{\mathrm{Raman}}$ yields $\Gamma_{\mathrm{Raman}}\approx 2\Gamma_{-+}^{\mathrm{Raman}}=5.0\times 10^{-4}T^{9}$. The prefactor in this expression corresponds to the value of the parameter $E_{t}=96$ introduced in Sec. V.1. The latter corresponds to $\bar{v}_{t}\equiv v_{t_{1}}v_{t_{2}}/\sqrt[5]{\left(v_{t_{1}}^{5}+v_{t_{2}}^{5}\right)/2}\approx 1.14\times 10^{3}\,\mathrm{m/s}$, a value close to the speed of sound extracted from the heat capacity data, $c=2\times 10^{3}$ m/s ($c$ is usually higher than the average transverse speed of sound $\bar{v}_{t}$). Using this value, at respectively 175 mT and 350 mT, our formula gives $\Gamma_{\mathrm{direct}}=0.9T$ and $14.9T$, which matches the order of magnitude of experimentally fitted dependencies, $\Gamma_{\mathrm{direct}}=7.2T$ and $47.1T$. In fact, this is a quite good agreement since the effect from other distortions is expected to contribute comparably to the rotational deformations. At 60 mT, our expression gives much smaller value, $\Gamma_{\mathrm{direct}}=0.013T$, than the experimentally fitted one $\Gamma_{\mathrm{direct}}=4.2T$. This discrepancy is easily understood given that at small values of applied field, the hyperfine/dipolar interactions and, accordingly, quantum tunneling effects influence considerably the total relaxation rate.

Interestingly, an improvement to the calculated results is achieved by a simple addition of a mean internal magnetic field resulting from nuclear spin and dipolar magnetic moments of surroundings molecules. In Fig. 4b, we present the computed total relaxation rate with an additional internal field of 115 mT. As can be seen, taking into account the internal field substantially boosts the overall agreement of the calculated results, $\Gamma_{\mathrm{direct}}=0.9T$, $7.0T$, and $46.3T$, with experimentally fitted $\Gamma_{\mathrm{direct}}=4.2T$, $7.2T$, and $47.1T$ at $H=60$ mT, 175 mT, and 350 mT, respectively. The worst agreement (but still close), as expected, is obtained for 60 mT. As already explained, this may result from the neglected effect of quantum tunneling as well as from the oversimplified mean field approximation, which does not account for the dynamics of the internal field. The latter, as well as the neglected effects from the distortions of the molecule, explains why such a high internal field, exceeding several times the expected one in molecules crystals Gatteschi et al. (2006), is needed to obtain the the quantitative agreement.

In order to check the field dependence of the relaxation rate predicted by the present theory, in Fig. 5 we show this dependence calculated with above expressions vs. experimental data for the same parameters. Quite similarly, without an additional internal magnetic field, the theoretical results at large fields yield a good agreement with the experiment for the orders of magnitude but become worse at small fields. However, with an addition of an internal field of 115 mT , the agreement between theory and experiment becomes striking with only a little difference in the small-field region as expected from the exclusion of the quantum tunneling effects and due to the approximation of a static internal field.

We note that a rationalization of these data was attempted earlier Gómez-Coca et al. (2014) within a mechanism involving the interaction of the electronic $S=3/2$ with the nuclear spin of Co(II). The authors of Ref. [Gómez-Coca et al., 2014] advocated in addition the hyperfine-nuclear interaction with an unrealistically high coupling constant ($\alpha_{I}=0.05\,\mathrm{cm^{-1}}$). The obtained agreement with experimental data (Fig. 5b in [Gómez-Coca et al., 2014]) is less perfect than in Fig. 5b, implying that the relaxation via acoustic phonons is probably more relevant for the presented compound.

VII Summary and discussions

In summary, we have derived the following new expressions for the low-temperature relaxation rate in the rotation approximation for the spin-phonon coupling:

• •

The transition rate between two states $\ket{\pm}$ of a ground Kramers doublet

• –

The direct process:

[TABLE]

• –

The Raman process:

[TABLE]

• •

The transition rate between two tunnel-split states $\ket{\pm}$ of a ground non-Kramers doublet

• –

The direct process:

[TABLE]

• –

The Raman process:

[TABLE]

In the case of Kramers doublets, all formulae are expressed via the $g$-factor of the ground doublet, $g_{\alpha},\alpha=x,y,z$, and the Zeeman splittings for the field applied in the corresponding direction, $W$ and $\Delta_{x,y}$. Eqs. (109) and (110) deliberately emphasize the $z$-direction, along the main $g$-factor, having in mind application to strongly axial magnetic centers ($g_{z}\gg g_{x},g_{y}$), otherwise being completely symmetric w.r.t $g_{\alpha}$ and $\Delta_{\alpha}$. In this way, the fourth-power dependence on the Zeeman splitting of the relaxation rate for a direct process Abragam and Bleaney (1970) appears in Eq. (109) in the form of three contributions involving the “bias” Zeeman splitting $W=g_{z}H_{z}$and the Zeeman “tunneling gap” $\Delta=\sqrt{\Delta_{x}^{2}+\Delta_{y}^{2}}.$ In the case of non-Kramers doublets, Eqs. (111) and (112), $W$ is the only Zeeman splitting allowed by the Griffith theorem Griffith (1971) and $\Delta$ is the intrinsic tunneling gap. Certainly, all results are independent of our choice of basis as physically expected.

These expressions have been derived under the assumption: $\Omega\ll\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}\ll\hbar\omega_{\xi 1}$ where $\xi$ denotes excited doublets. The remarkable feature of these expressions is their universal form for both Kramers and non-Kramers system. Furthermore, these universal formulae are entirely expressed via measurable or ab initio computable physical quantities, which facilitates their use.

It is worthy to note that although these results were derived supposing that the $g$-tensor of the multiplet is isotropic, an extension to a general $g$-tensor is trivial. In fact, the main difference now is in the matrix representation of the operator $S_{\alpha}$ (or $J_{\alpha}$), $\alpha=x,y,z$, in the sub-basis of the ground doublet. In particular, for Kramers system, this can be found by simply solving the system of equation $\mu_{\alpha}=-\frac{1}{2}g_{\alpha}\sigma_{\alpha}=-\sum_{\beta}g_{\alpha\beta}S_{\beta}$ where $g_{\alpha}$ are the principal values of the ground doublet’s $g$-tensor and $g_{\alpha\beta}$ are component of the multiplet’s $g$-tensor in the chosen reference frame. Notice that a new derivation also requires a change in the Zeeman Hamiltonian to $\mathcal{H}_{Z}=\sum_{\alpha,\beta}H_{\alpha}g_{\alpha\beta}S_{\beta}$ 555We consider multiplets $S$ or $J$ not far from isotropic case, when the Zeeman operator contains the first-rank contribution only Chibotaru (2013). This is by far the most common case.and accordingly the use of new expression $\left[\mathcal{H}_{Z},\bm{\mathrm{S}}\right]=i\left(\bm{\mathrm{S}}\times\bm{\mathrm{H}}^{\prime}\right)$ where $H_{\alpha}^{\prime}=\sum_{\beta}H_{\beta}g_{\beta\alpha}$. In practice, SMMs are usually characterized by a spin $S$ ($J$) with an almost isotropic $g$-tensor. Therefore, for the sake of simplicity we do not show here the results for this general case.

The total relaxation rate is then simply,

[TABLE]

Since the total relaxation rate was derived in the eigenstates basis, the phonon-induced quantum tunneling effect has been included from the start Ho and Chibotaru (2017). Hyperfine and dipolar interaction have been excluded from the present treatment, so the obtained results are strictly applicable to diluted samples.

The obtained expressions for direct process transition rate can be applied to tunnel-split states of any excited doublet at any temperature as well. This is because the derivation of the direct process transition rate for both Kramers and non-Kramers doublets was in fact not restricted by low temperature. However, the same cannot be said about the Raman process since its derivation was essentially done under the condition $\hbar\omega_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}\ll\hbar\omega_{\xi 1}\,\forall\ket{\xi}\neq\ket{\pm}$, which is typically satisfied at low temperature but generally not fulfilled at high temperature.

Within the rotational approximation to the spin-phonon coupling, one of the notable findings about the Raman transition rates is that the first-order process is always dominated by the second-order one for both Kramers and non-Kramers systems. Furthermore, whereas the Raman process in a Kramers system, as usual, is subject to the $T^{9}$ dependence emerging from the degeneracy of the doublet the Raman process in a non-Kramers system does not show the expected $T^{7}$ dependence but the $T^{9}$ one. Our derivation in fact demonstrates that this property results from the destructive interference of the main contributions to the $T^{7}$ dependence of the first-order and second-order Raman process, which are together left with the $T^{9}$ dependent terms only.

Moreover, contrary to the traditional belief that the Raman process in non-Kramers system is independent of the field Abragam and Bleaney (1970), we reveal that in fact, at least within the contribution from rotation spin-phonon coupling, the Raman process is strongly dependent on the applied field. Furthermore, it behaves in a similar manner as the incoherent quantum tunneling rate in the sense that it is significantly reduced when the system is out of resonance Gatteschi et al. (2006); Prokof’ev and Stamp (1998, 2000); Vijayaraghavan and Garg (2009); Garanin and Chudnovsky (1997); Leuenberger and Loss (2000). Not less interesting, for Kramers rigid molecules, the Raman process only depends on the magnetic field orientation but independent of its magnitude. Consequently, the Raman process keeps unchanged under a variation of the magnetic field magnitude. These properties, luckily, are easy verifiable in crystals of rigid SMMs. It should be noted that these results are only valid in the rotational approximation for the spin-phonon coupling. Although this approximation is sufficient to describe qualitatively the temperature and field dependence of the relaxation rate in the low-$T$ domain, to achieve a quantitative level of description, the contribution of intramolecular distortions, of the nuclear spins and dipolar interactions should be included. In some situation, when the latter become important, they can modify qualitatively the predictions of the present model. For instance, the field dependence of the relaxation rate for the Kramers doublets, predicted here as a monotonic function (Fig. 2c), does show a maximum in the Er-trensal complex Pedersen et al. (2014), due to an important contribution from the incoherent tunneling induced by nuclear spin and dipolar interactions

Regarding the question of the temperature dependence of the total spin-lattice relaxation rate, the present study shows that the spin-phonon interaction at low temperature always gives rise to a *$T$- *dependent relaxation. This is physically obvious given the fact that the phonon density varies strongly with energy. As a result, our finding rules out the common belief of a $T$-independent phonon-induced quantum tunneling rate at low temperature and zero field, and hence leaves the explanation of $T$-independence of the relaxation rate observed experimentally to other relaxation mechanisms, such as hyperfine Prokof’ev and Stamp (1996, 2000); Garanin et al. (2000) or dipolar interaction Prokof’ev and Stamp (1998); Fernández and Alonso (2003); Cuccoli et al. (1999).

A new strategy in designing SMMs with long relaxation times may also be drawn from our estimation on the direct and Raman process relaxation rates shown in Sec. V. Practically, it appears that the relaxation rate caused by the rotational deformations is quite small, even when the system’s anisotropy is not high in the case of Kramers system, or the intrinsic tunneling splitting is considerable in the case of non-Kramers system. This observation thus suggests that besides focusing on systems with strong anisotropy, it would be equally efficient to suppress the vibrational deformations by devising mechanically rigid SMMs, provided the quantum tunneling caused by hyperfine/dipolar interaction is suppressed.

In conclusion, we have derived universal “parameter-free” analytical expressions for the spin-phonon relaxation rate at low temperature for crystals of rigid molecules (both Kramers and non-Kramers) with arbitrary ZFS or crystal-field Hamiltonian, using a rotational approximation for spin-phonon interaction. These expressions offer a handy and quick estimation of the relaxation rate in arbitrary systems. Since the rotation deformation is always present in spin-phonon interaction, the relaxation rates found here set a lower bound on the total relaxation rates.

Acknowledgements.

L. T. A. H. is grateful to Dr. N. Iwahara and Dr. L. Ungur for helpful discussions and would like to acknowledge financial support from the Flemish Science Foundation (FWO).

Appendix A Averaging relaxation process rates over phonon modes

A.1 Direct process

We start by simplifying the expression (17) for the direct process transition rate $\Gamma_{-+}^{\mathrm{dr}}$ using the Debye dispersion relation and the fact that $\left[\bm{\mathrm{k}}\times\bm{\mathrm{e}}_{\bm{\mathrm{k}}t}\right]=\pm k\bm{\mathrm{e}}_{\bm{\mathrm{k}}t^{\prime}}$, where $t,t^{\prime}$ can be $t_{1}$ or $t_{2}$, to obtain

[TABLE]

Replacing $\sum_{\bm{\mathrm{k}}}$ by $\frac{V}{8\pi^{3}}\int d^{3}k$, then substituting $k=\omega_{\bm{\mathrm{k}}t}/v_{t}$, gives rise to

[TABLE]

The term inside the square bracket is the average of $\left(\bm{\mathrm{\Xi}}\cdot\bm{\mathrm{e}}_{\bm{\mathrm{k}}t}\right)\left(\bm{\mathrm{\Xi}}\cdot\bm{\mathrm{e}}_{\bm{\mathrm{k}}t}\right)^{*}$ over the unit $k-$sphere. Given the uniform distribution of the polarization vector $\bm{\mathrm{e}}_{\bm{\mathrm{k}}t_{1}}$ ($\bm{\mathrm{e}}_{\bm{\mathrm{k}}t_{2}}$) within the $k$-sphere, the integration over the spherical coordinates $\left(\theta_{\bm{\mathrm{k}}},\varphi_{\bm{\mathrm{k}}}\right)$ of $\bm{\mathrm{k}}$ can be changed to the spherical coordinates of vector $\bm{\mathrm{e}}_{\bm{\mathrm{k}}t_{1}}$ ($\bm{\mathrm{e}}_{\bm{\mathrm{k}}t_{2}}$) itself. The above integral then transforms into

[TABLE]

where $\bm{\mathrm{n}}$ is a normal unit vector on a sphere’s surface. Substituting Eq. (116) into Eq. (115) yield Eq. (18).

A.2 Raman process

Starting with the form (61) for $\Gamma_{-+}^{\mathrm{Raman}}$, valid for both Kramers and non-Kramers systems, averaging over all transversal phonons modes for this form can be done by first taking average over the phonon numbers $n_{\bm{\mathrm{k}}\lambda_{\bm{\mathrm{k}}}}$ and $n_{\bm{\mathrm{q}}\lambda_{\bm{\mathrm{q}}}}$, then substituting the Eq. (20) for $M_{\mathrm{ph}}^{\alpha\beta}$ into Eq. (61) gives:

[TABLE]

In analogy to the direct process, Eq. (115), using the Debye dispersion relation and the property of transverse phonons $\left[\bm{\mathrm{k}}\times\bm{\mathrm{e}}_{\bm{\mathrm{k}}t}\right]=\pm k\bm{\mathrm{e}}_{\bm{\mathrm{k}}t^{\prime}}$, we have

[TABLE]

which is Eq. (62).

Appendix B $\left(2S-1\right)\times\left(2S-1\right)$ submatrix diagonalization

Here we prove that a diagonalization of the submatrix of dimension $\left(2S-1\right)\times\left(2S-1\right)$ corresponding to the excited doublets $m^{\mathrm{th}}$ $\left(m\neq 1\right)$ can be avoided for the perturbed states $\ket{+}$ and $\ket{-}$. Since $\mathcal{H}_{Z}$ is small, a diagonalization of this matrix will lead to the following correction for $\ket{m}$,

[TABLE]

and similarly for $\ket{-_{m}}$, where $\ket{\pm_{m}}$ are eigenstates of the $m^{\mathrm{th}}$ doublet in the presence of the field. Since only the zeroth order term plays a role in the expression for $\ket{\pm}$ of the ground doublet, we can neglect the third term in the r.h.s. of Eq. (118). The corrections to the states $\ket{\pm}$ then become:

[TABLE]

The diagonalization of the submatrix block $\left(2\times 2\right)$ corresponding to the $m^{\mathrm{th}}$ doublet keeps the completeness of the basis. In other words, we have,

[TABLE]

$\ket{\pm}$ then become,

[TABLE]

Appendix C Time-reversal symmetry and matrix elements of time-even/time-odd operators

The derivation of the relaxation rate in the main text has made use of the following relations between matrix elements imposed by the time reversal symmetry:

• •

$\bm{\mathrm{\mathcal{H}}}^{\left(1\right)}\equiv i\left[\mathcal{H}_{A},\bm{\mathrm{S}}\right]$ is a time-even operator.

• •

Time-even operators matrix elements relation in a Kramers system: $\mathcal{H}_{\alpha\bar{\beta}}=-\mathcal{H}_{\beta\bar{\alpha}}$ , $\mathcal{H}_{\bar{\alpha}\beta}=-\mathcal{H}_{\bar{\beta}\alpha}$, $\mathcal{H}_{\bar{\alpha}\bar{\beta}}=\mathcal{H}_{\beta\alpha}$, and $\mathcal{H}_{\bar{\alpha}\alpha}=\mathcal{H}_{\alpha\bar{\alpha}}=0$.

• •

Time-odd operators matrix elements relation in a Kramers system: $O_{\alpha\bar{\beta}}=O_{\beta\bar{\alpha}}$, $O_{\bar{\alpha}\beta}=O_{\bar{\beta}\alpha}$, $O_{\bar{\alpha}\bar{\beta}}=-O_{\beta\alpha}$, and $O_{\bar{\alpha}\bar{\alpha}}=-O_{\alpha\alpha}$.

• •

Time-even operators matrix elements relation in a non-Kramers system: $\mathcal{H}_{\alpha\bar{\beta}}=\mathcal{H}_{\beta\bar{\alpha}}$, $\mathcal{H}_{\bar{\alpha}\beta}=\mathcal{H}_{\bar{\beta}\alpha}$, $\mathcal{H}_{\bar{\alpha}\bar{\beta}}=\mathcal{H}_{\beta\alpha}$, and $\mathcal{H}_{\bar{\alpha}\bar{\alpha}}=\mathcal{H}_{\alpha\alpha}$.

• •

Time-odd operators matrix elements relation in a non-Kramers system: $O_{\alpha\bar{\beta}}=-O_{\beta\bar{\alpha}}$, $O_{\bar{\alpha}\beta}=-O_{\bar{\beta}\alpha}$, $O_{\bar{\alpha}\bar{\beta}}=-O_{\beta\alpha}$, $O_{\bar{\alpha}\bar{\alpha}}=-O_{\alpha\alpha}$, and $O_{\alpha\bar{\alpha}}=O_{\bar{\alpha}\alpha}=0$.

These relations are easily proved by following the steps described in Chapter 15 of Abragam and Bleaney’s book Abragam and Bleaney (1970).