Spin relaxation: under the sun, anything new?*

Bogdan A. Rodin1,2,3 and Daniel Abergel1 1Laboratoire des biomolécules, LBM, Département de chimie, Ecole normale supérieure, PSL University, Sorbonne Université, CNRS, 75005 Paris, France 2International Tomography Center, Siberian Branch of the Russian Academy of Science, Novosibirsk 630090, Russia 3Novosibirsk State University, Novosibirsk 630090, Russia Correspondence: Daniel Abergel (daniel.abergel@ens.fr)


Introduction
Relaxation is the process through which a system loses energy to its environment to eventually reach a state of thermal equilibrium. Spin-lattice relaxation has been described so as to describe the way spins transfers energy to orientation degrees of freedom. Since the early days of NMR, it was the subject of numerous studies, both theoretical and encompassing a wide range of domains of applications. NMR relaxation theory has really been contemporary of the early days of NMR, and it was 10 formalized by several of its "founding fathers" (Bloembergen et al. (1948); Bloch (1957Bloch ( , 1956; Abragam (1961) ;Redfield (1957Redfield ( , 1965). This is a rather usual approach in the theory of open systems, which has been widely used in various domains of Physics (VanKampen (1981)). In this perspective, relaxation is the result of the dynamical coupling of a small ensemble of spins (the "system") coupled to a large ensemble of particles, or "degrees of freedom" (the "lattice") that is at thermal (Boltzmann) equilibrium and is endowed with an infinite heat capacity, thereby constituting a thermal reservoir. This very general 15 approach has led to many theoretical predictions with far reaching practical applications in the domain of magnetic resonance spectroscopy. In particular, the role played by molecular motions has been put in good use to extract dynamical information on complex molecular objects. Thus, with the development of "spin engineering" techniques, it has become possible to measure selected relaxation rates with high accuracy in a broad range of problems, with the prospect of relating such observables to models of molecular dynamics. 20 It has been shown recently by Bengs and Levitt (Bengs and Levitt (2020)) that the formulation of relaxation currently used by NMR spectroscopists may lead to erroneous predictions in the case of a two-spin system prepared in a singlet state. Such an unexpected behaviour was ascribed to the fact that some of the assumptions of the theory may not be fulfilled in NMR * This paper is dedicated to Prof. Geoffrey Bodenhausen on the occasion of his 70th birthday.
Taking the derivative (and renaming time indices), one gets: The dynamics restricted to the spin system is obtained by eliminating the bath variables. This is achieved by performing a partial trace over the bath degrees of freedom: 70 σ(t) = T race B {ρ(t)} = T race B e −iHTt ρ(0)e iHTt = T race B e LTt ρ(0) Hence, the spin density operator in the interaction representation: obeys: 75 Initially, the system and the bath are assumed to be completely decorrelated: where, the exact ρ e B denotes the bath density operator in thermal equilibrium. With these assumptions, one has: In the latter expression, the term T race B {L * 1 (t 1 )ρ * (0)} may be assumed to be zero or can be incorporated in the main system 80 liouvillian L S (Abragam (1961); Redfield (1957)). If the spin density operator is assumed to only moderately depart from its initial state, the zero-order approximation σ(t) ≈ σ(0) can be retained, so that σ(0) can be replaced by σ(t) in Eq. 11 (Abragam (1961); Redfield (1957)): This evolution equation can be transformed back to the Schrödinger representation. Since: we have: Thus, using the property: and using the definition e LSt σ * (t) = σ(t), one obtains the master equation in the Schrödinger representation: The term T race B {L * 1 (t)L * 1 (t )ρ e B } in equation 19 is a correlation operator acting on the spin system. It projects the system-bath (spin-lattice) coupling onto the bath degrees of freedom. This spin operator therefore carries the statistical properties of the bath, described by its equilibrium, stationary, density operator. Finally, assuming that the correlation operator decays to zero in a time τ c much shorter than the period over which the density matrix varies significantly, the upper limit of the integral can be extended to +∞ and: or, in the hamiltonian representation: One therefore obtains a master equation of the Redfield kind: ]} is the Redfield (relaxation) operator.

Formulation of the master equation in operator form
In spin relaxation theory, it is customary to express the relaxation equation in operator form, which often provides a clearer representation of the spin-bath coupling dynamics. Here, the coupling Hamiltonian is assumed to have the form of a sum of 110 terms, each of which factorizes into a product of lattice B q and spin S q operators.
with the interaction representation:

115
Using results of the preceding section, the Redfield equation writes: Each term in the sum writes:

120
where the notation: has been introduced. The B q (t)B q (t ) e are the bath (lattice) correlation functions, and in contrast to equation 20, these denote functions rather than operators.
Using the conventional decomposition of the spin operators into a sum of eigenoperators of the liouvillian L S = [H S , •]: one has: 130 one obtains from equation 29: Intoducing the secular approximation ω q p + ω q p = 0, so that p = p , q = −q , and renaming indices, this reduces to: The assumption that the bath is in a stationary state, [F, ρ e B , ] = 0, confers some properties to the correlation functions. Thus, the bath correlation functions are also stationary. Indeed, one has: In addition, because T r(AB) * = T r(B † A † ), it is easy to show that: Besides, using the property that B −q = B q † :

The Redfield equation is equivalent to the Lindblad form of the relaxation equation
It is now straightforward to show that the conventional Bloch-Redfield-Abragam perturbative approach of relaxation is equiv-150 alent to the Lindblad formulation of dissipative systems. Indeed, changing indices in the first term and using the property , and setting τ → −τ one has, from Eq. 33: Using the stationarity property (equation 34) of the bath correlation functions leads to: so that: where the "right" spectral density J q,−q R (ω q p ) of the bath is given by: {·, ·} denotes the anti-commutator. The operator appearing in the right handside of Eq. 40 is the Lindblad dissipation superoperator as: so that equation 40 writes: 170 (2007)), which is thus derived from the usual quantized theory of relaxation (Bloch (1957); Redfield (1957); Hubbard (1961); Abragam (1961)).

Equation 43 (as well as 40) is a Lindblad equation (Lindblad (1976); Alicki and Lendi
In principle, one should not expect the perturbative approach to yield an irreversible dissipative operator that is equivalent to a Lindblad operator. Indeed, this equivalence requires the Markovian and the short correlation time assumptions that make the evolution equation depend on the density operator at the present time t only, not on its previous history. Moreover, it also 175 requires the secular approximation that eliminates the time dependence of the spin operators of the coupling hamiltonian in Eq. 29, leading to Eq. 33. The combination of these conditions lead to the semi-group property and the Lindblad form of the relaxation operator.

An alternative formulation (equivalent to Lindblad)
It is possible to obtain an alternative and completely equivalent form of the Redfield equation. Expanding Eq. 33, one gets: where β = /kT . The latter equations express a Kubo kind of relation (Kubo (1957)): given in the Appendix (see Eqs A1 and A2). After some straightforward manipulations one obtains: Thus, collecting and rearranging terms, one gets: The Boltzmann factor can be expanded in series: The first term on the right hand side of this equation is the usual double commutator, while the second term represents the 195 thermal effect of the Boltzmann equilibrium of the lattice. This term, equal to 1−exp (−βω q p ), vanishes for infinite temperature.

A Pseudo-classical version of the Redfield equation
A pseudo-, or semi-classical version of the master equation can be extremely useful, allowing one to make use of models derived in the framework of classical mechanics to calculate spectral density functions. A general definition of the classical correlation function of two dynamical variables A and B is (Evans and Moriss (2008)): where the brackets indicate classical ensemble average.In the case of stationary processes, the following properties of a correlation functions can be deduced. Its complex conjugate C * AB (t) is therefore (Evans and Moriss (2008)): For an autocorrelation function of A, C AA (t), one has: Eq. 51 shows that, in the general case where the autocorrelation function is complex: ) and C i AA (t) = Im(C AA (t)) are even and odd functions of time, since C * . This implies that the associated spectral density, A semi-classical relaxation theory should provide spectral density functions obeying the general classical mechanics requirements detailed above. It is clear, however, that in the quantum case, where A and B are in general non commuting operators, the above symmetry relations do no apply. It is nevertheless possible to define a symmetrized correlation function as: which is real (when C AB (t) is stationary). Note also that the bath operator correlation functions have the property: (correlation functions are assumed stationary), so that: Using the definitions C −q,q , the relations Eqs. 54 and 55 show that the aver- ) obey the classical correlation function property C q * (τ ) = C q (−τ ).
A semi-classical version of the Redfield equation is thus obtained by using the spectral density function J q (ω) obtained from the Fourier transform of the symmetrized correlation function C q (τ ): Using Eqs. 56 and 46, it is therefore possible to derive an alternative expression of the master equation. This gives: or: Finally, collecting and rearranging terms, one gets: The final relaxation superoperator, which defines the relaxation of density matrix asσ * S (t) =Γσ * (t), may be written as: HereD β is the thermalized double commutator superoprator, that composes the superoperator from the the two operators A, B in the way: where f β (ω) = 2 1+e βω , and the dot is the place for operator on which superoperator is applied. It is easy to see that when Lindblad dissipator is still easily recognizable, as one may substitute 56 to 43 and get: Equation 59  The condition of Eq. 51 then implies that the spectral density function J q (ω) is even:

Simplifications in the high temperature approximation
When the largest eigenvalue of the operators S q p is such that max(βω q p ) 1, Eq. 59 takes the simpler form: where {., .} denotes the anticommutator. The "low order" approximation The assumption that the density operator is always close to the fully disordered state: ||σ − 1 A || 1, where A is the dimension of the density operator, was made by Redfield (Redfield (1957)). Limiting the expansion to zeroth order, Eq. 62 becomes (Hubbard (1961)): Moreover, using the property, Eq. 30, and the Taylor expansion of the exponential, it is straighforward to show that: Therefore, when the density operator is in thermal equilibrium determined by the hamiltonian H S , σ eq = T race(exp −βHS ) −1 exp −βHS , one can show that Rσ eq = 0, where R is defined by Eq. 59. Then, discarding terms that are second order or higher in max(βω q p ) 1, in Eq. 62, one gets the semi-classical formulation of the Redfield equation (Abragam (1961)):

270
The evolution of the expected value of an operator is given by the alternative master equation: In this expression, the second term on the right hand side contains the "thermal" contributions to relaxation, and can be selectively neglected for terms that are higher than first order in ω q p kT . That is, each term in the development such that: 275 at all times can be discarded. Eq. 67 is in principle a less stringent condition and may provide criteria for the quasi-, pseudoclassical approximation -a test that can be verified a posteriori. It may thus provide a way to select which parts of the density operator can be discarded (neglected) and which must be retained in order to get an approximate analytical solution.
The simple case of a two-spin system "Double commutator" versus "thermal" contributions 280 The differences of contributions between the first ("double commutator") and the second ("thermal") series of terms in equation 66 are illustrated in Figures 1 and 2 in the case of a pair of like spins 1 2 subject to relaxation caused by a mutual dipolar (dipole-dipole -DD) interaction and the presence of a randomly fluctuating field (ran). Simulations were performed assuming Lorentzian spectral density functions: with correlation times τ ran = 60 ps and τ C = 8 ps, for the random field and dipolar interactions, respectively. The dipolar coupling constant b 12 = −(µ 0 /4π)γ 2 I r −3 12 refers to the dipole-dipole coupling constant, γ I the gyromagnetic ratio and r 12 the internuclear distance. In these simulations, b 12 = 35 · 10 3 Hz. These values were chosen to give T 1 ≈ 2 s and T S ≈ 20 s. This is the case for both relaxation, dipolar and random field fluctuations, mechanisms (Figs 1(a) and (b)). Moreover, the values 295 chosen for the simulation imply that the dipolar contribution (of the double commutator in this case) to the total relaxation rate <Ȯ > (t) is much larger than than the one of the random field.
The situation is strikingly different when the spins are initially prepared in a singlet order. Here, the thermal correction (blue) is negligible with respect to the double commutator (red) contribution to the rate of change <Ȯ > (t) only for the dipole-dipole mechanism ( Fig. 1(c)). In contrast, for terms originating from random field relaxation, both thermal and dipolar terms are of 300 comparable orders of magnitude ( Fig. 1(c)). These are the terms that cannot be neglected in an approximate solution of Eq. 66. Fig. 1(c) also shows that the dipolar contribution to <Ȯ > (t) increases with time, which is consistent with the progressive depletion of the singlet order (immune to dipolar relaxation). And for random relaxation, which mainly affects the singlet order in this example, Fig. 1(d) illustrates the fact that the weight of the thermal contribution decays with time, with the concomitant increase of the double commutator term, which is also due to the progressive depletion of the singlet order.

305
The situation depicted in Fig. 2 is different, and shows the the rate of change of the expected value < O > (t) where O = 1 4 (I · S), for the same initial state conditions as above. In this case, the dipole-dipole simply does not contribute to <Ȯ > (t), as expected from symmetry considerations. This is of course the case whatever the initial state (inverted magnetization (Fig.2(a)) or singlet order (Fig.2 (c)). This illustrates the known fact that singlet state is immune to dipolar relaxation for symmetry reasons.

310
Alternatively, when the spins are prepared in the −(I z + S z ) state, both thermal and double commutators contribute, albeit to a negligible amount, showing that the spins evolve mostly towards magnetization (compare the scales with Fig. 1(b)), and that only a negligible part is transferred to singlet order.
Interestingly, Fig. 2(d) shows that there is no thermal contribution (blue) to the rate <Ȯ > (t), and that, starting from a singlet order, its evolution can be predicted by discarding the thermal thermal terms of Eq. 66, and therefore retaining the 315 simple, double commutator, expression for the relaxation master equation.

Singlet-triplet conversion
The recent achievement of the Lindblad approach was the description of the magnetization relaxation of a two-spin system prepared in a singlet state (Bengs and Levitt (2020)). In that paper, detailed balance was enforced through the Schofield  Table B1). (a) and (b) correspond respectively to the dipolar and random field relaxation of spins inverted from a Boltzmann equilibrium; (c) and (d) correspond respectively to the dipolar and random field relaxation of spins initially prepared in a singlet state.
procedure (Schofield (1960)), whereby spectral density functions are built from classical ones through the transformation: where J cl (ω) refers to the classical spectral density function. Eq. 69 is one among several that have been proposed to make classical spectral density functions asymmetric so as to obey the detailed balance condition (White et al. (1988); Egorov and Skinner (1998);Egorov et al. (1999); Ramirez and López-Ciudad (2004);Frommhold (1993)). In Ref. (Bengs and Levitt (2020)), J(ω) was assumed to be any kind of spectral density function obtained through classical models, such as diffusion 325 jumps, etc. The distinction between left and right spectral densities that appears in the course of the conventional perturbative derivation of the master equation, was not made there. Moreover, the detailed balance condition appears as an additional requirement, as this condition is not implied by Linblad's approach that merely provides the general mathematical structure of the evolution equation obeyed by the density operator that complies with the requirements of quantum mechanics, in the presence of a Markovian dissipative process (Lindblad (1976); Alicki and Lendi (2007)).

330
In the following, we derive the evolution of the magnetization of a two-spin system using the singlet-triplet population basis and compare results obtained by both approaches. As above (and in Ref. (Bengs and Levitt (2020)) the relaxation superoperator Γ is the sum of contributions from mutual dipole-dipole relaxation (Γ DD ) and the interaction with a partially correlated random field (Γ ran ): The irreducible tensor operator T λµ representation is better suited to derive analytical solutions for the problem at hand, wherê Γ is expressed, according to Eq. 61, as: , T 2µ ], The T λµ are eigenoperators of the main Zeeman Hamiltonian and their expressions are shown in Appendix B. ω (i) rms is the root mean square fluctuation of the random field acting on spin I i , and in the isotropic case considered here, is identical for both 340 nuclei, so that ω rms . The coefficient −1 ≤ κ 12 ≤ 1 describes the degree of correlation of random field fluctuations on the 1 and 2 nuclei. By definition, κ 11 = κ 22 = 1. In order to simplify the notations, we will henceforth drop the subscript (κ 12 → κ).
In the extreme narrowing regime where ωτ C,ran 1, the spectral densities become frequency independent and J(ω) ≈ 2τ .
The dipole-dipole and random field contributions to the longitudinal relaxation rate constant R 1 = R DD 1 + R ran 1 are given by, 345 according to Eq. 61: where f β (ω) = 2 1+e βω .It is interesting to note that in this model the relaxation rates do not depend on the temperature, in contrast to Ref. (Bengs and Levitt (2020)). This is not due to any approximation, rather, it arises from the fact that f β (ω) + 350 f β (−ω) = 2, which is explicit in equation 72. Similarly, the singlet order relaxation rate is given by: , which does not depend on the temperature. For sake of comparison with the results of (Bengs and Levitt (2020)), we use the singlet-triplet population basis, where the singlet and triplet states are defined as: where |α and |β denote the Zeeman spin states of an isolated spin-1/2 nucleus with z-projection of +1/2 and -1/2. In this representation, the population block of relaxation superoperator 71 is given by: where Σ i denotes the sum of the terms alongside respective column. This matrix is very similar to one introduced in (Bengs and Levitt (2020)), but with the substitution of term θ(ω) = exp (−βω/2) by f β (ω).
In the high-temperature limit, both terms become approximately equal, θ(ω) ≈ f β (ω) ≈ 1 − βω 2 and the difference between θ(ω) and f β (ω) becomes significant only in case of extremely low temperatures or very high frequencies. As could be expected 365 in this limit, the time evolution of the z magnetization is given by the same biexponential behavior as in Ref. (Bengs and Levitt (2020)): with the same coefficients A 1 = R S 2(R 1 − R s ) and A S = −2R 1 + R S 2(R 1 − R s ) . The fact that A 1 and A S are the ones found in (Bengs and Levitt (2020)) is expected, because in this limit R 1 and R S do not depend on the temperature factor. In the usual conditions 370 of high but not infinite temperature, it is found that the next nonzero terms in the expansions of A 1 (β) and A S (β) are of degree β 2 , and therefore do not contribute in the regime where ωβ 1. These expressions can be found in Appendix C.
The foregoing discussion has shown that in the high temperature approximation, the exact "thermalization" procedure of the spectral density function is irrelevant, as all models are equivalent in these conditions. Indeed, in a field of 23.5 T (1000 MHz resonance proton frequency) the temperature at which ω H ≈ kT , where both approaches may lead to significant differences, 375 is T ≈ 50 mK. These are unrealistic experimental conditions. Alternatively, the master equation of Eqs. (59-61) may well be of use in the context of DNP to describe the electron spin-lattice relaxation outside of the high-temperature limit, through direct spin-phonon coupling at temperatures below 4 K and at high fields, where this process is predominant.
6 Conclusion: a remark on the "semiclassical" theory In the semiclassical viewpoint (as in Ref. (Abragam (1961)), for instance), the effect of the bath is taken into account through a stochastic spin hamiltonian, the spatial part of which is a function of the lattice variables and is a random function of  (Bengs and Levitt (2020)) showed that the usual semi-classical master equation was not able to 385 predict the correct magnetization evolution of a two-spin system prepared in a singlet state. The authors' interpretation was that "thermalization" was not properly taken into account by the inhomogeneous semiclassical version of the Redfield equation, which was ascribed to the fact that the classical approach does not comply with the detailed balance condition. In fact, detailed balance is statistical by nature, and per se, compatible with a semiclassical approach.
In the semi-classical theory of NMR relaxation, it is assumed without justification that the correlation function is an even 390 function of the time (Abragam (1961); Redfield (1957)), which indeed implies that detailed balance is not obeyed, as in this case the spectral density function J(ω) is an even function of frequency. If, on the other hand, detailed balance is taken into account, so that J(−ω) = e − ω/kT J(ω) and the general relations Eq. 49 or 51 obeyed by correlation functions are retained, it is easy to show from the symmetry properties of the spectral density function that the semiclassical Redfield equation (Abragam (1961)): is obeyed, with this definition of J(ω). However, the expected equilibrium density operator is not a stationary state of Eq. 77 in this case, which illustrates the (also known) fact that this condition alone is insufficient to completely determine the transition probabilities of the bath in the absence of a dynamical model for the latter. On the other hand, it is possible to describe the dynamics of a classical system where microscopic irreversibility, i.e., detailed balance, is ensured. This is straightforward from 400 the definition of the correlation function of a phase variable in classical mechanics, A(t)B * = dqdpρ e Be iLt A, where L is the classical Liouvillian acting on the phase space (Evans and Moriss (2008)). In addition, general procedures have been used that provide Fokker-Planck or master equations for diffusion that obey the detailed balance condition, yielding classical spectral density functions that comply with the Boltzmann equilibrium distribution and the classical laws of motion of the bath (see for instance (VanKampen (1981); Risken (1972); Wassam et al. (1980))), in particular in the context of magnetic resonance 405 (Stillman and Freed (1980)). Note that in the fully quantum approach, detailed balance is ensured by assuming that the bath is in a stationary state defined by a Boltzmann distribution of its energy states. Thus, the "irreducible" difference between the semi-classical and the fully quantized theory lies in the fact that the bath operators do not mutually commute, which prevents the expression in Eq. 27 to reduce to the double commutator. Both thermodynamic and quantum mechanical effects are thus entangled in the fully quantum mechanical treatment of relaxation.
since B −q = B q † . Similarly, one has: Noting that | f |B −q |f | = |B −q f f | = |B q † f f | = |B q * f f | = |B q f f |, and exchanging indices f ↔ f , one gets, from Eq. A2: Appendix B: Eigenoperators for a homonuclear coupled spin 1/2 pair The case of homonuclear spin pair the main hamiltonian is defined as: where ω 0 = −γB, and γ is the magnetogyric ratio and B is the field strength. The eigenoperators for this hamiltonain are summarized in the where the coefficients A1 and AS was shown before in the main text, and the coefficients C1 and CS is defined as: