Effect of phonons on the electron spin resonance absorption spectrum

The unavoidable presence of vibrations in solid-state devices can drastically modify the expected electron spin resonance (ESR) absorption spectrum in magnetically active systems. In this work, we model the effect of phonons and temperature on the ESR signal in molecular systems with strong $E \otimes e$ Jahn-Teller (JT) effect and an electronic spin-$1/2$. Our microscopic model considers the linear JT interaction with a continuum of phonon modes, the spin-orbit coupling, the Zeeman effect, and the response of the system under a weak oscillating magnetic field. We derive a Lindblad master equation for the orbital and spin degrees of freedom, where one- and two-phonon processes are considered for the phonon-induced relaxation, and the thermal dependence of Ham reduction factors is calculated. We find that the suppression of ESR signals is due to phonon broadening but not based on the common assumption of orbital quenching. Our results can be applied to explain the experimentally observed absence of the ESR signal in color centers in diamond, such as the neutral nitrogen-vacancy and negatively charged silicon-vacancy color centers in diamond.


Introduction
Since its discovery in 1945 [1], electron paramagnetic resonance (EPR) or electron spin resonance (ESR) has been extensively used to characterize molecular systems [2], quantum dots [3], metal complexes [4], organic radicals [5] or defects in solid-state systems [6]. Commonly, ESR consists of applying a constant frequency microwave field while a magnetic field is swept across spin-flip resonance transitions. Crucial information has been obtained in several biochemical problems [7], on Hyperfine interactions with high resolution [8] and on magnetic ions in metals [9], to name a few examples. On the other hand, by analyzing the ESR absorption spectrum and dispersion lines, effects such as anisotropy of gyromagnetic factors [10], Ham reduction factors [11,12], and spin relaxation rates [13] have been studied.
In many systems were a paramagnetic behavior is expected, however, no ESR absorption is detected [14]. Especially on those systems that are strongly coupled to vibrations [15], the ESR response is not observed at high temperatures. Particularly, experimental measurements with defects in diamond reveal an absence of the ESR signal at room temperature, which is observed for the negatively charged silicon-vacancy (SiV − ) [16] and the neutral nitrogen-vacancy (NV 0 ) centers in diamond [14]. In both cases, spin S = 1/2 quantum systems dominated by the dynamical Jahn-Teller (JT) effect are quoted as examples of ESR suppression due to the fast orbital averaging effect (Ham reduction factors), but with no mention of a theoretical explanation. Here, we demonstrate that a strong dynamical JT effect primarily modifies the orbital states through the Ham reduction factors, which changes the energy of the resonant transitions of the system instead of the amplitude of the ESR signal. Only by including phononinduced relaxation processes we reproduce the observed suppression of the ESR signal at room temperature, which is conclusively explained by the phonon broadening effects induced by two-phonon Raman processes.
Although many works have been devoted to phenomenologically describing the ESR response, still few of them deal with microscopic aspects of the electron-phonon coupling and its thermal effects. Several authors have addressed the problem of modeling the ESR response using out-of-equilibrium models [17,18,19]. Historically, the linear response theory associated with irreversible processes was introduced in order to give a more satisfactory description of the magnetic resonance phenomenon [19]. Nevertheless, the inclusion of thermal effects induced by electron-phonon processes on the ESR response has not yet been microscopically solved in trigonal systems with a strong E ⊗ e JT effect. In this work, we present a microscopic model for the ESR response in a generic E ⊗ e ⊗ SU (2) system, where we include the spin-1/2 degree of freedom. Such a system can be, for instance, the NV 0 and the SiV − centers, or more general: the negatively charged group IV-vacancy centers. Furthermore, we derive a thermal dependence of Ham reduction factors, which is consistent with recent ab initio calculations for defects in diamond [33] (zero temperature) and the original theory presented by Ham [11,12]. This paper is organized as follows. In section 2, we introduce the Hamiltonian of the E ⊗e⊗SU (2) system, where a large band gap between the ground and first-optically excited state is assumed. In section 3, we derive an effective Hamiltonian for the orbital and spin degrees of freedom induced by the strong linear JT coupling. Furthermore, we introduce an analytical expression for the thermal dependence of Ham reduction factors, where the zero-temperature case is obtained from recent ab initio calculations.
In section 4, we introduce the phonon-induced relaxation processes using a Lindblad master equation that includes one-and two-phonon processes. In section 5, we use the linear response theory, and we find an expression for the ESR absorption spectrum in terms of the resonant frequencies and the phonon relaxation rates, where a super-Ohmic environment is assumed. Finally, we discuss the suppresion of the ESR signal, and in section 7, we summarize our results.

Hamiltonian of the system
Let us consider a trigonal system with D 3 symmetry whose degenerate orbital states |X and |Y with symmetry E are coupled to multiple e-phonon modes. Traditionally this type of orbital-phonon interaction is modeled using the E ⊗ e JT theory [20]. Now, we shall include an internal spin degree of freedom S = 1/2 to model the effect of phonons on the ESR signal. Towards this end, we incorporate the spin states |↑ = |m s = +1/2 and |↓ = |m s = −1/2 into the dynamics of the system. In general, under the presence of a constant B and weak time-dependent b(t) magnetic fields, the Hamiltonian of the E ⊗ e ⊗ SU (2) system is given bŷ For multiple phonon modes the E ⊗ e JT Hamiltonian is defined as [20] where E 0 is the energy of the degenerate states |X, ↑ , |X, ↓ , |Y, ↑ and |Y, ↓ . The last two terms of the right-hand of Hamiltonian (2) describe the linear electron-phonon coupling and the energy of the harmonic oscillators. Here, k ∈ E means that we are considering phonon modes belonging to the irreducible representations E, and in what follow we called e-phonon modes. The term F k is the linear vibronic coupling constant, while µ k and ω k are the reduced mass and the phonon frequencies, respectively. In the orbital manifold {|X , |Y }, the operators are defined asσ x = |X Y | + |Y X|, σ y = −i|X Y | + i|Y X| andσ z = |X X| − |Y Y |. The second term in equation (1) is the spin-orbit interaction Acoustic Phonons weak coupling Jahn-Teller Figure 1. Schematic representation of the E ⊗ e ⊗ SU(2) system perturbed by an oscillating magnetic field. The dashed box enclosed the usual E ⊗ e Jahn-Teller system, which is coupled to the spin states via the spin-orbit coupling constant λ. The static magnetic field interacts with both orbital and spin degrees of freedom, while the oscillating magnetic field only perturbs the spin states. Resonant acoustic phonons are depicted as a thermal bath interacting with the orbital states.
where λ > 0 is the spin-orbit coupling constant, L is the orbital angular momentum operator, and S = 1 2 (ŝ x ,ŝ y ,ŝ z ). The spin operators are given byŝ x = |↑ ↓| + |↓ ↑|, s y = −i |↑ ↓| + i |↓ ↑| andŝ z = |↑ ↑| − |↓ ↓|. In the term −λL zŜz we are assuming that orbital states with symmetry different from E are far in energy from the doubly degenerate E-states. As a consequence, we neglect the effect of operatorsL x andL y in our manifold since these orbital operators mix orbital states separated by a large band gap (ground and excited states). Additionally, we haveL z = σ y which is a useful representation for the z-component of the angular momentum in the orbital basis {|X , |Y }. The third term in equation (1) is the spin and orbital Zeeman Hamiltonian associated with the static magnetic field B = (B x , B y , B z ) ( = 1) where γ s = µ B g s / and γ L = µ B g L / are the spin and orbital gyromagnetic constants, respectively, with µ B being the Bohr magneton. Here, g s = 2 and g L = 1 are the Landé g-factors. Since we have neglected the operatorsL x andL y it follows that L · B ≈L z B z , which justifies the approximation on equation (4). The fourth term in equation (1) is the strain Hamiltonian for a trigonal system [20,21] where γ x and γ y are the coupling constants that describe the effect of static distortions. Finally, the last term on equation (1) is the time-dependent Hamiltonian associated with the weak oscillating magnetic field b(t) = b 1 cos(ωt)e x Figure 1 shows a schematic illustration of the E ⊗ e ⊗ SU(2) system, illustrating the dynamics presented in our model. In the next sections, we derive the effective dynamics for the orbital and spin degrees of freedom, and we include the effect of a phononic bath composed of acoustic phonons.

Effective dynamics for the orbital and spin degrees of freedom
In this section, the important case of color centers in diamond is analysed to emphasizes the role of a strong JT effect on the orbital and spin dynamics. First, we focus on the characteristic time scales of the JT effect in comparison with the phonon-induced relaxation processes, which will be our most important starting point. The JT effect has been numerically studied using first-principles calculations for the 3 E excited state of the NV − center [22], both excited and ground states of the SiV 0 [30] and SiV − centers [33]. We focus our analysis to the particular case of the SiV − center in diamond to illustrate the main ideas of the microscopic model. For the SIV, the numerically estimated JT and barrier energies are E JT = 42.3 meV and δ JT = 3.0 meV, respectively [33]. Using the previous values and the relations E JT = F 2 /(2/( ω − 2G)) and δ JT = 4E JT G/( ω + 2G) it follows that the linear vibronic JT coupling is F = 83.34 meV, which is associated to a local vibrational mode with energy ω = 85.2 meV. The latter implies that the time scale associated with the JT interaction is τ JT ∼ g −1 JT ≈ 0.3 ps, which is a fast dynamics. We expect a similar time scale for the JT dynamics (few picoseconds) for the ground state of the neutral nitrogen-vacancy center. On the other hand, experimental measurements of the phonon-induced relaxation rates show that Γ ph ∼ (10 −5 − 10 7 ) s −1 for temperature ranging from few mK to room temperature [23,24,25,26,27,28], where a microscopic model for the electron-phonon processes of the SiV − center is discussed in Ref. [29]. As a consequence, phonon-induced relaxation processes are relevant for times larger than τ ph−relax ∼ 0.1 µs. Therefore, in any physical system satisfying the condition τ JT τ ph−relax the strong JT effect modifies the energy levels of the system in the short-time dynamics, and the phonon-induced relaxation processes will be resonant to an effective Hamiltonian dressed by the JT effect, which is the original idea of the Ham reduction factors [11,12].
In the time scale 0 < t ∼ τ JT , and assuming that τ JT τ ph−relax , the system is ruled by the following Liouvillevon Neumann equation whereρ t is the density matrix associated with the orbital, spin, and e-phonons degrees of freedom. Here, the time-dependent magnetic field is not considered as it has a much slower dynamics. In equation (7), we have an approximate Hamiltonian because we are neglecting the electron-phonon interaction responsible for the phonon-induced relaxation processes. Now, we will focus our analysis on the case of strong JT interaction. In such a case, the HamiltonianĤ JT is not a small perturbation, and therefore the usual weak-coupling approximation of the master equation cannot be applied to derive the dynamics for the orbital and spin states. An alternative approach is to find a new reference frame or unitary transformation in which the coupling with e-phonons will be small, and then we could apply perturbation theory. To accomplish the previous observation, we introduce the following unitary transformation whereT k =P kxσz +P kyσx is a linear phonon-orbital operator and k = −F k /µ k ω k , where F k , µ k and ω k are the parameters introduced in the JT Hamiltonian (2). The summation introduced in equation (8) only considers the contribution of e-phonon modes responsible of the JT effect. The transformation (8) was originally introduced by Ham in his pioneering work of the quenching effect of the orbital operators in E ⊗ e systems with strong JT interaction [12]. The coefficient k ∼ F k / ω k can be understood as the ratio between the vibronic coupling and the phonon energy, and for color centers in diamond this factor is small ( ≈ 0.1 for the SiV − center). In the regime F k /µ k ω k 1, we obtainĤ The first two terms of the right-hand side describe the free energy of the degenerate orbital states and e-phonon modes, respectively. More importantly, the last term defines an exchange of angular momentum between orbital states and e-phonons. The orbital and phonon momentum operators are defined asL z = σ y andL ph kz = −1 (P kxQky − P kyQkx ), respectively. After applying the transformation (8) the Liouvillevon Neumann equation (7) To leading order in k and calculatingρ s = Tr ph (ρ t ), we obtain the following effective closed dynamics for the orbital and spin degrees of freedom The most important result is the effective Hamiltonian where p and q are the Ham reduction factors which satisfy the relation q = (1 + p)/2 [12]. In figure 2 we illustrate the energy levels of the SiV − center as a function Figure 2. Energy levels of the E ⊗ e ⊗ SU (2) system as a function of the magnitude of the magnetic field B for a fixed temperature. The filled coloured areas represent the possible values of the energies E i for 0 ≤ θ ≤ π/2, λ = 50 GHz and γ x,y = 2 GHz. The transitions between the states are illustrated with black arrows which can be realized by one-and two-phonon processes.
of the magnetic field B = B(cos φ sin θe x + sin φ sin θe y + cos θe z ), where the energy gap p √ λ 2 + Υ 2 is affected by the Ham reduction factor p, with Υ = γ 2 x + γ 2 y being the total strain coupling. In the weak coupling regime ( k 1) and to first-order, the Ham reduction factor is reduced to where g k = F k /2µω k , k B is the Boltzmann constant and T is the reservoir temperature. For a detailed derivation of the Ham reduction factor see Appendix A.

Thermal dependence of the Ham reduction factor
To calculate the thermal dependence of the Ham reduction factor p for a continuum of ephonon modes is necessary to consider corrections beyond the first-order approximation given in equation (12). In this direction, second-order corrections lead to [31] From equations (12) and (13) we observe that the Taylor expansion of the Ham reduction factor has the form p . Therefore, when e-phonons are considered as a continuum field, we obtain the following analytical expression where J(ω) is the phononic spectral density function associated with e-phonon modes In principle, the spectral density function takes into account the contribution of all e-phonon modes associated with the JT effect and must be numerically calculated. Recent ab initio calculations show no indication of the presence of quasi-localized phonon modes with symmetry e g in the ground state of the SiV − center [32]. Assuming the absence of distinct quasi-localized resonances in the shape of the phononic spectral density function J(ω), we expect a dominant contribution coming from lattice phonons. Therefore, as a first theoretical approach, we consider the following super-Ohmic model for the phononic spectral density function where α = Ωg 2 JT /(2π 2 v 3 s ω D ) and ω c 2 THz is the cut-off frequency for the acoustic phonons in diamond. For a diamond lattice we have Ω = a 3 , where a = 3.57Å is the lattice constant, v s = 1.2×10 4 m/s is the speed of sound, ω D ≈ 40×10 12 Hz is the Debye frequency, and g JT is the electron-phonon coupling constant. By direct integration of equation (14) using the spectral density function J acous (ω) (16), we obtain where . Numerical calculations of the Ham reduction factor at zero temperature for the SiV − center in diamond leads to p num ≈ 0.308 [33]. By comparison we deduce that g JT /2π = (v 3 s 2 ω D ln(p −1 num )/(4Ωω 2 c )) 1/2 which can be used to calculate the electron-phonon coupling associated with e-phonons for a fixed set of parameters (Ω, ω c , v s , p num ).
The cut-off frequency ω c is a lattice-dependent parameter, and therefore it is instructive to model the expected thermal dependence of p acous (T ) for systems having a different ω c but equal initial value p(0) ≈ 0.3. In figure 3 we show the thermal dependence of the Ham reduction factor for different situations that fulfil the condition p acous (0) = 0.3. We observe that the contribution of multiple e-phonon modes can drastically decrease the Ham reduction factor. This is one of the most important results of this work, and we will study further implications of this behavior on the ESR signal. It is important to observe what happens in the high-temperature regime and weak static strain, i.e. γ x,y λ. This is the case for color centers in diamond with a small strain induced by the crystal distortion. At high temperatures (T ∼ 300 K) the Ham reduction factor p in equation (11) eliminates the contribution of the spin-orbit interaction leading to an effective HamiltonianĤ eff ≈ E 0 + γ s S · B. Thus, under these approximations, in an ESR experiment at room temperature, the resonances will be centered around the Zeeman frequencies. However, under realistic experimental conditions, the peaks of the ESR signal are broadened by the phonon relaxation rates. In the next section, we introduce the effect of phonon relaxation in order to quantify the role of one-and two-phonon processes on the ESR absorption spectrum.

Open dynamics and phonon relaxation processes
To include the phonon-induced relaxation processes, we consider the dynamics in the time scale t τ JT . For color centers in diamond, the contribution of one-and twophonon processes are crucial to model the thermal dependence of the phonon relaxation rates and their effects on the ESR signal. From experiments and theory is known that the most relevant phonon relaxation processes observed in color centers are the following: i) direct processes with the energy condition ω k = E i − E j , and ii) two-phonon processes [23,29,34]. In our model, the energy condition of the previous processes fully matches the energy gap E i − E j between different eigenstates of the effective Hamiltonian (11), where ω k are the acoustic phonon frequencies of the diamond lattice. As a consequence of the strong JT effect, the relaxation processes start to be observable in a time scale where the orbital reduction (Ham reduction factors) has been reached. Formally, the relaxation rates can be derived from the electron-phonon interaction and employing the Fermi golden rule theory [26]. In the basis spanned by the eigenstates of the effective Hamiltonian (11), we can write the electron-phonon Hamiltonian for one-and two-phonon processes as follows: where phonons with a linear dispersion relation ω k = v s |k| are considered in the electron-phonon Hamiltonian, being v s the speed of sound and k the wavevector for phonons. Here, λ ij,k and λ ij,kk are the linear and quadratic electron-phonon coupling constants, respectively andb † k (b k ) is the creation (annihilation) operator for acoustic phonons. The energy of the acoustic phonons is described by the Hamiltonian In order to derive the open dynamics, we consider acoustic phonons at thermal equilibrium and weakly coupled to the system. The latter is not a strong assumption and is usually satisfied for color centers in diamond for temperatures ranging from 10 mK to 500 K [26]. After applying the Born and secular approximations, we obtain the following Markovian master equation for the orbital and spin degrees of freedom whereL ij = |i j| defines a set of normalized projection operators, |i are the normalized eigenstates of the effective Hamiltonian (11), and γ ij are the phonon relaxation rates associated with the transition |i → |j . In the basis |i , the off-diagonal matrix elements can be analytically solved, and are given by ρ ij (t) = ρ ij (0)exp [− (iω ij + Γ ij ) t] for i = j, where ω ij = (E i − E j )/ and Γ ij is the effective relaxation rate To evaluate the individual rates we use the decomposition γ ij = γ 1−ph ij +γ 2−ph ij , where γ 1−ph ij and γ 2−ph ij are the relaxation rates induced by one-and two-phonon processes, respectively. Now, we discuss the low and high temperature regimes for the particular case of color centers. At low temperatures, T < 10 K, acoustic phonons dominate the thermal dependence of γ ij , and therefore γ ij ≈ γ 1−ph ij . In such a case, acoustic phonons can be described by a smooth phononic spectral density function, which scales as ∼ ω 3 in a three-dimensional lattice. By employing the Fermi golden rule theory to first-order using the linear contribution of the electron-phonon interaction (18), we obtain the following expression for the one-phonon relaxation rates where Ω is the volume of the unit cell in a diamond lattice, E = (F acous /ω ph ) 2 /(2µ), F acous is the coupling constant associated with acoustic phonons, and n(ω) = [exp( ω/k B T ) − 1] −1 is the mean number of phonons at thermal equilibrium. The expression for γ 1−ph ij has the same structure as the photon-induced relaxation rates of the Wigner-Weisskopf theory of the spontaneous emission. In our case, the rate Γ ij (21) will contribute to both absorption and emission processes, leading to a linear thermal scaling at high temperatures, i.e., Γ ij ∝ T [29] within our temperature range of analysis.
At high temperatures, T > 100 K, experimental measurements of the longitudinal relaxation time T 1 with negatively charged nitrogen-vacancy [23] and neutral siliconvacancy [34] centers shows that two-phonon processes are composed by two important contributions: i) Raman processes leading to 1/T 1 ∝ T n (n = 5 for NV − and n = 7 for the SiV 0 ) and ii) Orbach-type processes leading to 1/T 1 ∝ [exp(∆E/k B T ) − 1] −1 , where ∆E is the energy of a phonon (∆E = 73 meV for NV − and ∆E = 22 meV for SiV 0 ). However, recent experimental measurements for the SiV − show that in the temperature region between 500 mK and 2.3 K there is a more involved temperature dependence for 1/T 1 , which is beyond of the scope of this work. At room temperatures, the Raman processes are the most dominant contribution [23,34], which allows us to neglect the Orbach rates in a phenomenological model at high temperatures.
We can address the problem of modelling the contribution of two-phonon processes by considering a rate described by γ 2−ph ij = A ij T n , where A ij is a fit parameter and n = 5, 7 for the NV − , SiV 0 center, respectively. A detailed derivation of the term A ij T n is presented in Appendix C. Also, we note that the observed T 3 law for the phononinduced Raman processes of the SiV − [29] is only valid for the electronic transition between excited and ground states, as there only phonon modes with energy higher than the spin-orbit coupling contribute to the Raman transition. This case is crucially different from the internal ground dynamics presented in this work. We remark that a full derivation of the Orbach rates is beyond the scope of this paper, and we do not consider it since we are interested in the phenomenology thermal response of the ESR signal. In the next section, we introduce the analytical expression for the ESR absorption spectrum derived from the linear response theory.

ESR absorption spectrum
The response of the system when the oscillating magnetic field is swept across spin-flip resonance transitions depends on the intensity of the interaction HamiltonianV (t) = γ s b 1Ŝx cos(ωt). Nevertheless, in the weak driven limit, i.e., when |γ s b 1 | max|Ĥ eff | (whereĤ eff is the effective Hamiltonian (11)), the absorption spectrum can be defined as the imaginary part of the dynamical susceptibility [35,36] whereŜ x (0) is the initial x-component of the spin operator and the expectation value Ŝ x (t)Ŝ x (0) ss = Tr(Ŝ x (t)Ŝ x (0)ρ ss ) is calculated using the steady stateρ ss obtained from the quantum master equation (20). It is convenient to write the spin operatorŜ x (0) in the basis spanned by the eigenstates of the effective Hamiltonian (11), obtaininĝ where g ij = i|Ŝ x |j . The time-dependent spin operatorŜ x (t) is calculated using the following quantum dynamical map whereρ (i,j) s (t) is the solution of the master equation (20) starting from the initial conditionρ s (0) = |i j|. By calculating the expectation value Ŝ x (t)Ŝ x (0) ss using the steady stateρ ss = exp(−βĤ eff )/Tr(exp(−βĤ eff )), we get the following analytical expression for the ESR absorption spectrum The ESR absorption spectrum is a sum of Lorentzian functions with peaks at the resonant frequencies ω ij = (E i −E j )/ obtained from the effective Hamiltonian (11) and broadened by the phonon relaxation rates Γ ij given in equation (21). The Boltzmann distribution p i (T ) is given by In the next section, we study the combined effect of the Ham reduction factor p acous (T ) and the phonon relaxation rates Γ ij on the ESR absorption spectrum.

Motional and dynamical suppression effects on the ESR absorption spectrum
In order to better illustrate the combined thermal effects induced by the Ham reduction factor p acous (T ) and the phonon relaxation rates Γ ij on the ESR absorption spectrum we analyse the individual transition between two eigenstates of the effective Hamiltonian (11). Let us consider the transition |1 ↔ |4 , where |1 and |4 are the lowest and highest energy levels of the effective Hamiltonian, respectively. For the particular case of color centers in diamond, the spin-orbit coupling constant is larger than the strain parameters, and therefore we can assume that λ (γ 2 x + γ 2 y ) 1/2 . In such a case, the resonant frequency is given by ω 41 ≈ p acous (T )λ + γ s B z , where B = B zẑ is the static magnetic field. At high temperatures we have p acous (T ) ≈ 0 (see Figure 3) leading to ω 41 ≈ γ s B z , which is a drastic reduction of the resonant frequency induced by the quenched spin-orbit interaction. In figure 4-(a) we plot the frequency ω 41 = p acous (T ) + γ s B z as a function of temperature. At zero temperature, we have ω 41 ≈ 50 GHz, however when the temperature increases the resonant frequency is fully determined by the Zeeman energy term γ s B z . Therefore, the peak of the ESR signal around ω 41 moves as a function of temperature as shown in figure 4-(a). This motional effect of the ESR signal is a consequence of the strong JT interaction, which is encapsulated in the Ham reduction factor.
In figure 4-(b) we plot the ESR signal for the transition |1 ↔ |4 in terms of the detuning δ = ω − ω 41 and considering B = (0, 0, 1000) G, p acous (0)λ = 50 GHz and γ x = γ y = 1 GHz, and a phenomenological magnetic noise Γ mag = 3 MHz. By looking the shape of the ESR absorption spectrum we confirm that the intensity is drastically suppressed at high temperatures. This effect can be explained in terms of the phonon relaxation rate Γ 14 and the Boltzmann distribution factor p(T ) introduced in equation (26). From equation (26), we deduce that the intensity and the full width at half maximum (FWHM) of the ESR signal are given by Γ and 2|g 14 | 2 p 1 (T )/Γ, respectively at the resonant frequency ω = ω 41 . The phonon relaxation rate is given by Γ = Γ mag + Γ 14 , where the constant magnetic noise Γ mag describes the effect of magnetic impurities in the lattice and Γ 14 is calculated from equation (21). Because of the dominant contribution of two-phonon processes at high temperatures, we obtain a temperature scaling T 7 for the relaxation rate Γ 14 , which is shown in figure 4-(d) for different two-phonon coupling constants λ 00 ≈ λ 00ij (see Appendix C for further details). As a consequence, the strong thermal dependence Γ 41 ∝ T 7 can explain the reduction of the ESR signal when temperature increases and its absence at high enough temperatures.
In addition, we plot the intensity of the ESR signal at the resonant frequency ω = ω 41 in figure 4-(c). We numerically confirm that the intensity is reduced to zero above temperatures T > 100 K. In other words, the suppression of the ESR signal at room temperature is because of the activation of two-phonon processes and not for the effect of the Ham reduction factors induced by the strong JT effect. Towards this  and for different temperatures. We use a static magnetic field with components B x = B y = 0 and B z = 1000 G, a spin-orbit coupling p acous (0)λ = 50 GHz, and a local strain γ x = γ y = 1 GHz. (c) Maximum value of the ERS signal in comparison with its value at zero temperature I max /I max (T = 0). (d) Effective relaxation rate Γ 14 in logarithmic scale as a function of temperature and for different two-phonon coupling constants λ 00 ≈ λ 00ij , see Appendix Appendix C. The coupling constant λ 00 is taken equal to 2.4 µeV in order to obtain the same temperature dependence of the neutral silicon-vacancy center [34]. direction, recent experimental observations of the ESR transitions associated with the SiV − centers showed that the signal is observed at 100 mK [39] and also at 4 K [40] but not at room temperature [16]. Our model predicts that at low temperatures, two-phonon processes are not activated, resulting in an ESR signal with a high intensity, as shown in figure 4-(b). We remark that other similar systems may not present these motional and suppression effects at high temperatures if the set of parameters λ, γ x,y , Γ ij and p acous (T ) are different and/or the assumptions presented here are not satisfied for them. For instance, the vanadium in hexagonal silicon-carbide, a defect system with trigonal symmetry, spin-1/2, and ground state with dynamic JT effect still shows evidence for an ESR signal above 50 K [41,42,43].

Conclusion
In summary, we have presented a microscopic model to study the effect of phonons and temperature on the ESR absorption spectrum in a generic E ⊗ e Jahn-Teller system with residual spin S = 1/2. We derived the effective dynamics for the orbital and spin degrees of freedom by finding an analytical expression for the thermal dependence of the Ham reduction factor p. Furthermore, we included one-and two-phonon processes into the relaxation dynamics to analyse thermal effects on the ESR signal. Using our model and the SiV − parameters, we presented a physical explanation for the absence of contrast in the ESR response at high temperatures, and also we show that the signal is recovered at temperatures of the order of tens-hundreds mK. Interestingly, the latter is experimentally confirmed by recent experiments with color centers in diamond [39,37]. Most importantly, we demonstrated that a strong dynamical JT effect modifies the orbital operators in a short time scale, leading to an effective Hamiltonian whose energy levels are dominated by the Zeeman interaction at high temperatures. On the contrary, phonon-induced relaxation processes sharply decrease the amplitude of the ESR signal due to phonon broadening effects. The model can be used to characterize new spin-1/2 systems at low temperatures for metrology and quantum information applications, where a strong JT effect is present.

Appendix A. Ham reduction factors
In order to introduce the transformation given in equation (8) we first describe the effect of two independent modesQ x andQ y and after that we extend our results to a continuum of e-phonon modes. Firstly, let us consider the following unitary transformation where H sJT is the JT Hamiltonian for a two vibrational modes,T =P xσz +P yσx and = −F/µ ω. This allows us to decouple orbital and momentum operators, and the transformed JT Hamiltonian at first order in (weak coupling regime) will be [12] However, we must also apply this transformation to the other terms of the Hamiltonian (2), in particular, we have to apply it to terms that depend on orbital operatorsσ i (i = x, y, z). After averaging over the phonon bath we obtain where the subindex e represents an averaging over the electronic part of the states and p and q are the Ham reduction factors [12] which satisfy the relation q = 1 2 (p + 1). We get where P 2 = P 2 x + P 2 y . Using the trigonometric identity sin 2 x = (1 − cos 2x)/2, we obtain Introducing the quantization of phonons P x,y = i µ ω/2(b † x,y − b x,y ) and using the Fock representation |k = |k x , k y , we deduce that Now, if we calculate the expectation value using a thermal phonon state described by ρ ph = e −β /Z ph , where Z ph = Tr ph (e −βH ph ) and β = 1/k B T , we obtain and e −βH ph |k = e −2βn(ω) |k . Finally, the following Ham reduction factor is obtained where 2n(ω) + 1 = coth ( ω/2k B T ). To second order in the above expression is reduced to The generalization of p (2) to a continuum of phonon modes can be made by adding the contribution of all phonon modes, following this we recover our result given in equation (12).

Appendix B. Effective Hamiltonian
As previously we discussed in Appendix A we first consider the effect of two independent phonon modes. Let us consider the Liouville equatioṅρ whereĤ =Ĥ JT +Ĥ so +Ĥ z is the total Hamiltonian of the system, wherê is the JT Hamiltonian. We use the same unitary transformation introduced in equation (8) where = −F/µ ω. Now, we consider the Born or Hartree approximationρ(t) = ρ s (t) ⊗ρ ph , where phonons are considered in thermal equilibrium. Expanding in powers of , we obtainρ where we have used Tr(T ρ ph ) = 0 sinceT is a linear operator in P x and P y over the phonon sub-space. By taking the trace over the phonon degrees of freedom, we deducėρ where Tr ph ([Ĥ ,ρ ]) = Tr ph (Ĥ ρ ) − Tr ph (ρ Ĥ ) = Tr ph (Ĥ ρ ) − H.C.. By simplicity we calculate the first term contribution associated with each Hamiltonian Tr ph (Ĥ soρ ) = pĤ soρs , (B.8) Tr ph (Ĥ zρ ) = pĤ zoρs +Ĥ zsρs (B.9) where p is the Ham reduction factor,Ĥ zo is the orbital Zeeman interaction andĤ zs is the spin Zeeman interaction. Finally, we obtaiṅρ whereĤ eff = pĤ so + pĤ zo +Ĥ zs is the effective Hamiltonian.

Appendix C. Two-phonon processes
To model two-phonon processes, we use the electron-phonon Hamiltonian given in equation (18). In a three-dimensional lattice without cubic symmetry, we can model the coupling constants in the continuum limit as [38]: where ω D is the Debye frequency in diamond and λ 0ij , λ 00ij are coupling constants. To introduce the transition rate between eigenstates of the effective Hamiltonian we will focus on phonon processes satisfying the energy condition ω k − ω k = (E i − E j )/ , where E i are the eigenenergies of the effective Hamiltonian (11). The transition rate between two eigenstates can be defined as [26] γ i→j = k: ω k =vs|k| k : ω k =vs|k | Γ j,n k −1,n k +1 i,n k ,n k . (C.3) To second order in perturbation theory the Fermi golden rule states that Γ j,n l ,n l i,n k ,n k = 2π 2 V j,n l ,n l i,n k ,n k + m p,p V m,np,n p j,n l ,n l V i,n k ,n k m,np,n p E i,n k ,n k − E m,np,n p 2 × δ(ω ji + n l ω l + n l ω l − n k ω k − n k ω k ), (C. 4) where V j,n l ,n l i,n k ,n k = i, n k , n k |Ĥ e−ph |j, n l , n l , |n k represents a Fock state for the k-th vibrational mode and E i,n k ,n k = (ω i + n k ω k + n k ω k ), with E i = ω i . In the continuum limit the density of states (DOS) due to acoustic phonons is given by (C.5) where dΩ = sin θdθdφ is the differential solid angle, ω D = v s k D , D 0 = Ω/(2π 2 v 3 s ) and Θ(x) is the Heaviside function. Using the DOS and the intermediate phonon states |n p , n p = {|n k − 1, n k , |n k , n k + 1 }, we obtain the following transition rate γ i→j = a ij T 5 + b ij T 6 + c ij T 7 .
(C.6) Each coefficient is defined as where n(x) = [exp(x)−1] −1 and x ij = ω ij /(k B T ). By considering that the coupling with two phonons is described by a real parameter λ ijkk ∈ R, we conclude that b ij = 0. In such case, we have γ i→j = a ij T 5 + c ij T 7 . For a single SiV − center under the presence of a static magnetic field B = B zẑ , the allowed transitions induced by phonons are given by |1 ↔ |3 and |2 ↔ |4 . Therefore, the only non-zero linear electron-phonon parameters are λ 042 , λ 024 , λ 031 and λ 013 . As a consequence, we obtain a SiV i→j = 0 due to the fact that m λ 0jm λ 0mi = 0 for every i, j = 1, 2, 3, 4. Now, we assume that λ 00ij ≈ λ 00 . Second, due to the Ham reduction factors the terms x ij = ω ij /(k B T ) ≈ 0 for T > 100 K and low magnetic fields. Therefore, we obtain the following contribution to the effective transition rate Γ ij (see equation (21)) Γ ij ≈ 10380πD 2 0 |λ 00 | 2 k 7 B 7 ( ω D ) 2 × T 7 , (C.12) where we have numerically solved the integral I = x D 0 f (x) dx for f (x) = x 6 n(x)[n(x) + 1], obtaining a constant value I = 432.487 in the temperature interval (0 − 300) K. In order to compare our model to real systems, we extract the experimental value 1/T 1 = A Raman T 7 for the neutral silicon-vacancy center, where A Ramam = 5.0 × 10 −13 s −1 K −7 [34]. By assuming that λ 00ij ≈ λ 00 , we estimate λ 00 to be approximately 2.4 µeV.