Spin relaxation: is there anything new under the Sun?

Spin relaxation has been at the core of many studies since the early days of nuclear magnetic resonance (NMR) and the underlying theory worked out by its founding fathers. This Bloch–Redfield–Abraham relaxation theory has been recently reinvestigated (Bengs and Levitt, 2020) in the perspective of Lindblad theory of quantum Markovian master equations in order to account for situations where the widely used semi-classical relaxation theory breaks down. In this article, we review the conventional approach of quantum mechanical theory of NMR relaxation and show that, under the usual assumptions, it is equivalent to the Lindblad formulation. We also comment on the debate on semi-classical versus quantum versions of spectral density functions involved in relaxation.

1 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 as the way spins transfer energy to orientation degrees of freedom. Since the early days of nuclear magnetic resonance (NMR), it has been the subject of numerous studies, both theoretical and encompassing a wide range of domains of applications. NMR relaxation theory really was contemporary in the early days of NMR, and it was formalized by several of its founding fathers (Bloembergen et al., 1948; Bloch, 1957, 1956; Abragam, 1961; Redfield, 1957, 1965). This is a rather usual approach to the theory of open systems, which has been widely used in various domains of physics (VanKampen, 1981). From 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 are at thermal (Boltzmann) equilibrium and are endowed with an infinite heat capacity, thereby constituting a thermal reservoir. This very general 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 to 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.

It has been shown recently by Bengs and Levitt (2020) that the formulation of Redfield's semi-classical theory of relaxation widely used by NMR spectroscopists may lead to erroneous predictions in the case of a two-spin system prepared in high-order states, such as singlet spin states.

Such an unexpected behavior was ascribed to the fact that some of the assumptions of the theory may not be fulfilled in NMR systems, and the authors solved the problem by making use of the Lindblad operator, which is commonly used in the theory of open quantum systems to account for dissipative Markovian phenomena, i.e., relaxation processes. Among other properties, the structure of the Lindblad operator ensures that the fundamental properties of the density operator ρ (ρ, which is hermitian and definite positive, and Tr(ρ)=1) are preserved (Lindblad, 1976; Alicki and Lendi, 2007). This article has sparked renewed interest regarding the general theory of relaxation, as first elaborated by Bloch, and its connection with the Lindblad theory of quantum dissipative systems (Barbara, 2021). This approach has been generalized in several recent works (Nathan and Rudner, 2020; Maimbourg et al., 2021).

The traditional description of NMR relaxation relies on the description of both the spin system and the lattice as quantum systems, an approach that leads to the celebrated Redfield equation (Bloch, 1956; Redfield, 1965; Hubbard, 1961; see also Goldman, 2001, for a more recent account on spin lattice relaxation). As far as the description of the lattice is concerned, this approach is challenging, and actually untractable, as a quantized description of the degrees of freedom involved in molecular motions (multiple bond rotations, overall tumbling of a molecule, etc.) is not manageable in practice, even for small molecular systems. For this reason, an alternative relaxation theory, where the spins are treated as quantum objects and the lattice is described with classical functions of the lattice degrees of freedom, was developed. This semi-classical approach has far-reaching practical consequences, as spin relaxation can then, in principle, be described using classical models of dynamics for molecular motions, and it has been extensively used over the years to describe and interpret spin relaxation experiments. However, this semi-classical theory predicts a non-Boltzmann equilibrium density operator, which requires an ad hoc thermal correction to the relaxation equations (Redfield, 1965; Abragam, 1961). More elaborate attempts have been made to overcome this limitation by modifying the relaxation operator itself and to enforce a Boltzmann equilibrium of the density operator (Jeener, 1982; Levitt and Bari, 1992). However, in the traditional semi-classical NMR relaxation approach and its later modifications, the spin–lattice interactions are accounted for by a stochastic, fluctuating Hamiltonian. This has important consequences. First, in the conventional Abraham–Redfield approach, the statistical properties of the (classical) lattice are not constrained by a fluctuation–dissipation kind of theorem that would enforce a Boltzmann equilibrium distribution and, therefore, a lattice temperature. Such a constraint is therefore not applied on the spins. Second, a major difference between the quantum and classical theories of spin relaxation is rooted in the non-commutation of quantum bath operators of the spin–lattice coupling Hamiltonian, which confers particular properties to the spin correlation and spectral density functions that are absent in the semi-classical theory. Both aspects, quantum and statistical, are entangled, and the relations between quantum and classical correlation functions will therefore be discussed.

It is the purpose of this paper to re-investigate these old questions in order to trace the roles and the consequences of the various assumptions of the traditional approach to relaxation developed in the early days of NMR.

2 Theory of spin relaxation in a thermal bath: a short review

2.1 Derivation of the master equation

The derivation follows the lines of Abragam (1961), Bloch (1956), Redfield (1957), and Hubbard (1961). Consider a spin system in its environment. The Hamiltonian of the spin system H_S accounts for the interaction of the spins with the magnetic fields (Zeeman interaction and interaction with a radio frequency (rf) field) and non-dissipative spin–spin interactions (scalar J coupling and dipolar coupling). The dynamics of the lattice (the bath) are described by H_B, and the spin–bath interaction Hamiltonian is H₁, as follows:

$$\begin{array}{}\text{(1)}& H_{\text{T}}=H_{\text{S}}+H_{\text{B}}+H_{\mathrm{1}},\end{array}$$

and the dynamics of the {spin–bath} system are described by the Liouville equation as follows:

$$\begin{array}{}\text{(2)}& \dot{\mathit{\rho }}\left(t\right)=-i[H_{\text{S}}+H_{\text{B}}+H_{\mathrm{1}},\mathit{\rho }(t\left)\right].\end{array}$$

Alternatively, it can be written as follows:

$$\begin{array}{}\text{(3)}& \dot{\mathit{\rho }}\left(t\right)={\mathcal{L}}_{\text{T}}\mathit{\rho }\left(t\right),\end{array}$$

where the operator ${\mathcal{L}}_{\text{T}}=-i[H_{\text{T}},\cdot ]$ is the total Liouvillian. Introducing the density operator in the interaction representation of the density operator ${\mathit{\rho }}^{\ast }\left(t\right)=e^{-{\mathcal{L}}_{\mathrm{0}}t}\mathit{\rho }\left(t\right)$, where ${\mathcal{L}}_{\mathrm{0}}={\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}}$ is the unperturbed Liouvillian, one has the following:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\rho }}^{\ast }\left(t\right)=-i\left[H_{\mathrm{1}}^{\ast }\right(t),{\mathit{\rho }}^{\ast }(t\left)\right]={\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathit{\rho }}^{\ast }\left(t\right).\end{array}$$

The Liouville equation is integrated to the second order as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{\rho }}^{\ast }\left(t\right)& ={\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}_{\mathrm{1}}{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t_{\mathrm{1}}\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)\\ & +\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}_{\mathrm{1}}\underset{\mathrm{0}}{\overset{t_{\mathrm{1}}}{\int }}\text{d}{t}_{\mathrm{2}}{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t_{\mathrm{1}}\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t_{\mathrm{2}}\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)+\mathrm{\dots }\end{array}\end{array}$$

Taking the derivative, one obtains the following:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\rho }}^{\ast }\left(t\right)={\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right).\end{array}$$

Finally, by making the change in the variables $\mathit{\tau }=t-t^{\prime }$, one obtains the master equation for the total density operator as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\rho }}^{\ast }\left(t\right)={\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}\mathit{\tau }{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }(t-\mathit{\tau }){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right).\end{array}$$

The dynamics restricted to the spin system are obtained by eliminating the bath variables. This is achieved by performing a partial trace over the bath degrees of freedom as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathit{\sigma }\left(t\right)& ={\text{tr}}_{\text{B}}\left\{\mathit{\rho }\left(t\right)\right\}={\text{tr}}_{\text{B}}\left\{e^{-i{H}_{\text{T}}t}\mathit{\rho }\left(\mathrm{0}\right)e^{i{H}_{\text{T}}t}\right\}\\ & ={\text{tr}}_{\text{B}}\left\{e^{{\mathcal{L}}_{\text{T}}t}\mathit{\rho }\left(\mathrm{0}\right)\right\}.\end{array}\end{array}$$

Hence, from Eq. (6), the spin density operator in the interaction representation is as follows:

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}^{\ast }\left(t\right)={\text{tr}}_{\text{B}}{\mathit{\rho }}^{\ast }\left(t\right),\end{array}$$

which obeys the following:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\sigma }}^{\ast }\left(t\right)& ={\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t_{\mathrm{1}}\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)\right\}\\ & +\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)\right\}.\end{array}\end{array}$$

Initially, the system and the bath are assumed to be completely decorrelated, as follows:

$$\begin{array}{}\text{(11)}& \mathit{\rho }\left(\mathrm{0}\right)={\mathit{\rho }}_{\text{B}}^e\otimes \mathit{\sigma }\left(\mathrm{0}\right),\end{array}$$

where the exact ${\mathit{\rho }}_{\text{B}}^e$ denotes the bath density operator in thermal equilibrium. With these assumptions, one has the following:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\sigma }}^{\ast }\left(t\right)& ={\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t_{\mathrm{1}}\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)\right\}\\ & +\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{\mathit{\sigma }}^{\ast }\left(\mathrm{0}\right).\end{array}\end{array}$$

In the latter expression, the term ${\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t_{\mathrm{1}}\right){\mathit{\rho }}^{\ast }\left(\mathrm{0}\right)\right\}$ may be assumed to be zero or can be incorporated in the main system Liouvillian ℒ_S as follows (Abragam, 1961; Redfield, 1957):

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\sigma }}^{\ast }\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{\mathit{\sigma }}^{\ast }\left(\mathrm{0}\right).\end{array}$$

If the spin density operator is assumed to only moderately depart from its initial state, then the following occurs:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathit{\sigma }\left(t\right)-\mathit{\sigma }\left(\mathrm{0}\right)}{\mathit{\sigma }\left(\mathrm{0}\right)}}\ll \mathrm{1},\end{array}$$

where the density operator varies only slightly from its initial state so that σ^∗(0) can be replaced by σ^∗(t) in Eq. (12), as follows (Abragam, 1961; Redfield, 1957):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\sigma }}^{\ast }\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{\mathit{\sigma }}^{\ast }\left(t\right).\end{array}$$

The master equation in the Schrödinger representation can be obtained from Eq. (15) as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& {\displaystyle \frac{\text{d}}{\text{d}t}}{\mathit{\sigma }}^{\ast }\left(t\right)\\ & ={\displaystyle \frac{\text{d}}{\text{d}t}}{e}^{-{\mathcal{L}}_{\text{S}}t}\mathit{\sigma }\left(t\right)=-{\mathcal{L}}_{\text{S}}{e}^{-{\mathcal{L}}_{\text{S}}t}\mathit{\sigma }\left(t\right)+e^{-{\mathcal{L}}_{\text{S}}t}{\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)\\ & =\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{\mathit{\sigma }}^{\ast }\left(t\right).\end{array}\end{array}$$

Therefore, one has the following:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)& ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)\\ & +e^{{\mathcal{L}}_{\text{S}}t}\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{\mathit{\sigma }}^{\ast }\left(t\right).\end{array}\end{array}$$

Reverting to the Schrödinger representation ${\mathit{\sigma }}^{\ast }\left(t\right)=e^{-t{L}_{\text{S}}}\mathit{\sigma }\left(t\right)$, one obtains the following master equation in the Schrödinger representation:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}{\mathcal{L}}_{\mathrm{1}}^{\ast }(t^{\prime }-t){\mathit{\rho }}_{\text{B}}^e\right\}\mathit{\sigma }\left(t\right).\end{array}$$

A derivation is given in Appendix A for reference. The term ${\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}{\mathcal{L}}_{\mathrm{1}}^{\ast }(t^{\prime }-t){\mathit{\rho }}_{\text{B}}^e\right\}$ in Eq. (18) 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 and 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 +∞. As above, making the change in the variables $\mathit{\tau }=t-t^{\prime }$, one obtains the following:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}{\mathcal{L}}_{\mathrm{1}}^{\ast }(-\mathit{\tau }){\mathit{\rho }}_{\text{B}}^e\right\}\mathit{\sigma }\left(t\right),\end{array}$$

or, in the following Hamiltonian representation, the following:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)& =-i[H_{\text{S}},\mathit{\sigma }(t\left)\right]\\ & -\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{\text{tr}}_{\text{B}}\left\{[H_{\mathrm{1}},[H_{\mathrm{1}}^{\ast }(-\mathit{\tau }),{\mathit{\rho }}_{\text{B}}^e\mathit{\sigma }\left(t\right)\left]\right]\right\}.\end{array}\end{array}$$

One, therefore, obtains a master equation of the Redfield kind, as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)=-i[H_{\text{S}},\mathit{\sigma }(t\left)\right]+\mathcal{R}\mathit{\sigma }\left(t\right),\end{array}$$

where $\mathcal{R}•=-{\int }_{\mathrm{0}}^{+\mathrm{\infty }}\text{d}\mathit{\tau }{\text{tr}}_{\text{B}}\left\{[H_{\mathrm{1}},[H_{\mathrm{1}}^{\ast }(-\mathit{\tau }),{\mathit{\rho }}_{\text{B}}^e•\left]\right]\right\}$ is the Redfield (relaxation) operator.

2.2 Formulation of the master equation in operator form

In spin relaxation theory, it is customary to express the relaxation equation in the 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 terms, each of which factorizes into a product of lattice B^q and spin S^q operators.

$$\begin{array}{}\text{(22)}& H_{\mathrm{1}}=\sum _q{S}^q{B}^q,\end{array}$$

with the interaction representation as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }{B}^q\left(t\right)=e^{i{H}_{\text{B}}t}{B}^q{e}^{-i{H}_{\text{B}}t}\text{(24)}& {\displaystyle }{S}^q\left(t\right)=e^{i{H}_{\text{S}}t}{S}^q{e}^{-i{H}_{\text{S}}t}.\end{array}$$

Using results of the preceding section (Eqs. 16–18), the Redfield equation becomes the following:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _{q,q^{\prime }}{\text{tr}}_{\text{B}}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}{t}^{\prime }\left[S^q\right(t\left)B^q\right(t),\\ & \left[S^{q^{\prime }}\right(t^{\prime }\left)B^{q^{\prime }}\right(t^{\prime }),{\mathit{\rho }}_B^e{\mathit{\sigma }}_{\text{S}}^{\ast }(t\left)\right]].\end{array}\end{array}$$

Each term in the sum becomes the following:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& {\text{tr}}_{\text{B}}\left[S^q\left(t\right)B^q\left(t\right),\left[S^{q^{\prime }}\left(t^{\prime }\right)B^{q^{\prime }}\left(t^{\prime }\right),{\mathit{\rho }}_B^e{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\right]\\ & =\left[S^q\left(t\right),S^{q^{\prime }}\left(t^{\prime }\right){\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\langle B^q\left(t\right)B^{q^{\prime }}\left(t^{\prime }\right){\rangle }^e\\ & +\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S^{q^{\prime }}\left(t^{\prime }\right),S^q\left(t\right)\right]\langle B^{q^{\prime }}\left(t^{\prime }\right)B^q\left(t\right){\rangle }^e,\end{array}\end{array}$$

where the notation, as in the following:

$$\begin{array}{}\text{(27)}& \langle B^q\left(t\right)B^{q^{\prime }}\left(t^{\prime }\right){\rangle }^e={\text{tr}}_{\text{B}}\left\{B^q\left(t\right)B^{q^{\prime }}\left(t^{\prime }\right){\mathit{\rho }}_B^e\right\},\end{array}$$

has been introduced. The $\langle B^q\left(t\right)B^{q^{\prime }}\left(t^{\prime }\right){\rangle }^e$ are the bath (lattice) correlation functions, and in contrast to Eq. (19), these denote usual time correlation functions rather than operators.

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _{q,q^{\prime }}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}{t}^{\prime }\left[S^q\right(t),S^{q^{\prime }}(t^{\prime }\left){\mathit{\sigma }}_{\text{S}}^{\ast }\right(t\left)\right]\langle B^q\left(t\right)B^{q^{\prime }}\left(t^{\prime }\right){\rangle }^e\\ & -\sum _{q,q^{\prime }}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}{t}^{\prime }\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\right(t\left)S^{q^{\prime }}\right(t^{\prime }),S^q(t\left)\right]\langle B^{q^{\prime }}\left(t^{\prime }\right)B^q\left(t\right){\rangle }^e.\end{array}\end{array}$$

Using the conventional decomposition of the spin operators into a sum of eigenoperators of the Liouvillian ${\mathcal{L}}_{\text{S}}=[H_{\text{S}},•]$ as follows:

$$\begin{array}{}\text{(29)}& [H_{\text{S}},S_p^q]={\mathit{\omega }}_p^q{S}_p^q,\end{array}$$

one has the following:

$$\begin{array}{}\text{(30)}& S^q\left(t\right)=e^{i{H}_{\text{S}}t}{S}^q{e}^{-i{H}_{\text{S}}t}=\sum _p{S}_p^q{e}^{i{\mathit{\omega }}_p^{q}t},\end{array}$$

which one obtains from Eq. (28), with the change in the integration variable $\mathit{\tau }=t-t^{\prime }$ as follows:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _{q,q^{\prime },p,p^{\prime }}{e}^{i({\mathit{\omega }}_p^q+{\mathit{\omega }}_{p^{\prime }}^{q^{\prime }})t}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[S_p^q,S_{p^{\prime }}^{q^{\prime }}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\\ & \langle B^q\left(t\right)B^{q^{\prime }}(t-\mathit{\tau }){\rangle }^e{e}^{-i{\mathit{\omega }}_{p^{\prime }}^{q^{\prime }}\mathit{\tau }}\\ & -\sum _{q,q^{\prime },p,p^{\prime }}{e}^{i({\mathit{\omega }}_p^q+{\mathit{\omega }}_{p^{\prime }}^{q^{\prime }})t}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_{p^{\prime }}^{q^{\prime }},S_p^q\right]\\ & \langle B^{q^{\prime }}(t-\mathit{\tau })B^q\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_{p^{\prime }}^{q^{\prime }}\mathit{\tau }}.\end{array}\end{array}$$

Introducing the secular approximation ${\mathit{\omega }}_p^q+{\mathit{\omega }}_{p^{\prime }}^{q^{\prime }}=\mathrm{0}$, so that $p=p^{\prime },q=-q^{\prime }$, and renaming the indices reduces to the following:

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _{p,q}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[S_p^{-q},S_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\\ & \langle B^{-q}\left(t\right)B^q(t-\mathit{\tau }){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & -\sum _{p,q}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q,S_p^{-q}\right]\\ & \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}.\end{array}\end{array}$$

The assumption that the bath is in a stationary state, $[H_{\text{B}},{\mathit{\rho }}_B^e,]=\mathrm{0}$, confers some properties to the correlation functions. Thus, the bath correlation functions are also stationary. Indeed, one has the following:

$$\begin{array}{}\text{(33)}& \begin{array}{rl}& \langle B^q\left(t\right)B^{-q}(t+\mathit{\tau }){\rangle }^e\\ & ={\text{tr}}_{\text{B}}\left\{e^{i{H}_{\text{B}}t}{B}^q{e}^{-i{H}_{\text{B}}t}{e}^{i{H}_{\text{B}}(t+\mathit{\tau })}{B}^{-q}{e}^{-i{H}_{\text{B}}(t+\mathit{\tau })}{\mathit{\rho }}_B^e\right\}\\ & ={\text{tr}}_{\text{B}}\left\{e^{i{H}_{\text{B}}(t-\mathit{\tau })}{B}^q{e}^{-i{H}_{\text{B}}(t-\mathit{\tau })}{e}^{i{H}_{\text{B}}t}{B}^{-q}{e}^{-i{H}_{\text{B}}t}{\mathit{\rho }}_B^e\right\}\\ & =\langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e.\end{array}\end{array}$$

In addition, because $\text{tr}(AB{)}^{\ast }=\text{tr}(B^{†}{A}^{†})$, it is easy to show the following:

$$\begin{array}{}\text{(34)}& \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^{e*}=\langle B^q\left(t\right)B^{-q}(t-\mathit{\tau }){\rangle }^e.\end{array}$$

Besides, using the property that $B^{-q}=B^{q†}$, as in the following:

$$\begin{array}{}\text{(35)}& \begin{array}{rl}& \langle B^{-q}\left(t\right)B^q(t-\mathit{\tau }){\rangle }^e\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}\langle f\left|e^{i{H}_{\text{B}}t}{B}^{-q}{e}^{-i{H}_{\text{B}}t}\right|f^{\prime }\rangle \\ & \langle f^{\prime }\left|e^{i{H}_{\text{B}}(t-\mathit{\tau })}{B}^q{e}^{-i{H}_{\text{B}}(t-\mathit{\tau })}{e}^{-\mathit{\beta }{H}_{\text{B}}}\right|f\rangle \\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}\langle f\left|B^{-q}\right|f^{\prime }\rangle e^{ift}{e}^{-i{f}^{\prime }t}\langle f^{\prime }\left|B^q\right|f\rangle \\ & e^{i{f}^{\prime }(t-\mathit{\tau })}{e}^{-if(t-\mathit{\tau })}{e}^{-\mathit{\beta }f}\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}\langle f\left|B^{-q}\right|f^{\prime }\rangle \langle f^{\prime }\left|B^q\right|f\rangle e^{-i(f^{\prime }-f)\mathit{\tau }}{e}^{-\mathit{\beta }f}\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f^{\prime }|B^q|f\rangle {|}^{\mathrm{2}}{e}^{-i(f^{\prime }-f)\mathit{\tau }}{e}^{-\mathit{\beta }f},\end{array}\end{array}$$

where the |f〉 are the eigenstates of the bath Hamiltonian, H_B, and $\mathit{\beta }=\frac{\mathrm{\hslash }}{kT}$, with the convention ℏ=1. In these equation, the notation $L={\text{tr}}_{\text{B}}{e}^{\mathit{\beta }{H}_{\text{B}}}$ was introduced.

3 The Redfield equation is equivalent to the Lindblad form of the relaxation equation

It is now straightforward to show that the conventional Bloch–Redfield–Abraham perturbative approach of relaxation is equivalent to the Lindblad formulation of dissipative systems. Indeed, by changing indices in the first term, using the property ${\mathit{\omega }}_p^{-q}=-{\mathit{\omega }}_p^q$, and setting $\mathit{\tau }\to -\mathit{\tau }$, one has, from Eq. (32), the following:

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _{p,q}\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\text{d}\mathit{\tau }\left[S_p^q,S_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\\ & \langle B^q\left(t\right)B^{-q}(t+\mathit{\tau }){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & -\sum _{p,q}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q,S_p^{-q}\right]\\ & \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}.\end{array}\end{array}$$

Using the stationarity property (Eq. 33) of the bath correlation functions leads to following:

$$\begin{array}{}\text{(37)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _{p,q}\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\text{d}\mathit{\tau }\left[S_p^q,S_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\\ & \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & -\sum _{p,q}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q,S_p^{-q}\right]\\ & \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\end{array}\text{(38)}& {\displaystyle }\begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& \approx -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{p,q}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[S_p^q,S_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\\ & \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{p,q}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q,S_p^{-q}\right]\\ & \langle B^q(t-\mathit{\tau })B^{-q}\left(t\right){\rangle }^e{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }},\end{array}\end{array}$$

so that, in the following:

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& \approx -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{p,q}\left[S_p^q,S_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]J_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{p,q}\left[{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q,S_p^{-q}\right]J_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & =\sum _{p,q}{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\left(S_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\right(t)S_p^q\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{S_p^q{S}_p^{-q},{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right\}),\end{array}\end{array}$$

where the right spectral density $J_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)$ of the bath is given by the following:

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)& =\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\langle B^q{B}^{-q}\left(\mathit{\tau }\right)\rangle e^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & =\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{C}_{\text{R}}^{q,-q}\left(\mathit{\tau }\right)e^{-i{\mathit{\omega }}_p^q\mathit{\tau }},\end{array}\end{array}$$

where $\mathit{\{}\cdot ,\cdot \mathit{\}}$ denotes the anti-commutator. The operator appearing on the right-hand side of Eq. (39) is the Lindblad dissipation super-operator, as follows:

$$\begin{array}{}\text{(41)}& \widehat{\mathcal{L}\mathcal{D}}[A,B]=A•B-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(BA•+•BA\right),\end{array}$$

so that Eq. (39) becomes the following:

$$\begin{array}{}\text{(42)}& {\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)=\sum _{p,q}{J}_R^{q,-q}\left({\mathit{\omega }}_p^q\right)\widehat{\mathcal{L}\mathcal{D}}\left[S_p^{-q},S_p^q\right]{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right).\end{array}$$

Equation (42) (and Eq. 39) is a Lindblad equation (Lindblad, 1976; Alicki and Lendi, 2007), which is thus derived from the usual quantized theory of relaxation (Bloch, 1957; Redfield, 1957; Hubbard, 1961; Abragam, 1961).

The fact that this derivation leads to the Lindblad equation is not obvious. In principle, one should not expect the perturbative approach to yield an irreversible dissipative operator equivalent to a Lindblad operator. In fact, 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 and not in its previous history. Moreover, it also requires the secular approximation that eliminates the time dependence of the spin operators of the coupling Hamiltonian in Eq. (28), leading to Eq. (32). The combination of these conditions lead to the semi-group property and the Lindblad form of the relaxation operator.

Another point is worth mentioning. The properties of the correlation functions that emerge through this procedure reflect the properties of the bath operators of the spin–bath coupling Hamiltonian and, therefore, convey additional properties that are not implied by the structure of the Lindblad in Eqs. (39) and (42). Such properties, arising from the non-commutation of the bath operators and the fact that the lattice is always in a stationary Boltzmann equilibrium (detailed in Sect. 3.1 below and in Appendix B), ensure a detailed balance and that the stationary state of the spin system is compatible with the Boltzmann equilibrium of the lattice, that is, the lattice temperature. In other words, the perturbative approach leads to a master equation of the Lindblad form, with additional physical properties that bear constraints from the lattice.

3.1 An alternative formulation (equivalent to Lindblad)

It is possible to obtain an alternative and completely equivalent form of the Redfield equation. Expanding Eq. (32), one obtains the following:

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _q\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{S}_p^{-q}{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\langle B^{-q}\left(\mathit{\tau }\right)B^q\rangle e^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & +\sum _q\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{S}_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q\langle B^q{B}^{-q}\left(\mathit{\tau }\right)\rangle e^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & +\sum _q\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^{-q}\langle B^{-q}\left(\mathit{\tau }\right)B^q\rangle e^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & -\sum _q\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{S}_p^{-q}\langle B^q{B}^{-q}\left(\mathit{\tau }\right)\rangle e^{-i{\mathit{\omega }}_p^q\mathit{\tau }}.\end{array}\end{array}$$

As above, the left-hand side spectral density function is defined as follows:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{J}_{\text{L}}^{-q,q}\left(\mathit{\omega }\right)& =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\langle B^{-q}\left(\mathit{\tau }\right)B^q\rangle e^{-i\mathit{\omega }\mathit{\tau }}\\ & =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{C}_{\text{L}}^{-q,q}\left(\mathit{\tau }\right)e^{-i\mathit{\omega }\mathit{\tau }}=e^{-\mathit{\beta }\mathit{\omega }}{J}_{\text{R}}^{q,-q}\left(\mathit{\omega }\right),\end{array}\end{array}$$

where $\mathit{\beta }=\mathrm{\hslash }/kT$ (ℏ=1). The latter equations express a Kubo kind of relation (Kubo, 1957), where $J_{\text{L}}^{-q,q}\left(\mathit{\omega }\right)=e^{-\mathit{\beta }\mathit{\omega }}{J}_{\text{R}}^{q,-q}\left(\mathit{\omega }\right)$. A proof thereof is given in the Appendix (see Eqs. B1 and B3). After some straightforward manipulations, one obtains the following:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{S}_p^{-q}{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)e^{-\mathit{\beta }{\mathit{\omega }}_p^q}{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{S}_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^{-q}{e}^{-\mathit{\beta }{\mathit{\omega }}_p^q}{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{S}_p^{-q}{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right).\end{array}\end{array}$$

Thus, by collecting and rearranging terms, one obtains the following:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\left(\left[S_p^{-q},{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q\right]\right.\\ & \left.-\left[S_p^{-q},S_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]e^{-\mathit{\beta }{\mathit{\omega }}_p^q}\right).\end{array}\end{array}$$

The Boltzmann factor can be expanded in a series, as follows:

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left[S_p^q,{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\right]\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & \left[S_p^{-q},S_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}}{n\mathrm{!}}}(-\mathit{\beta }{\mathit{\omega }}_p^q{)}^n\right].\end{array}\end{array}$$

The first term on the right-hand side of this equation is the usual double commutator and the symmetry, while the second term represents the thermal effect of the Boltzmann equilibrium of the lattice. This term, equal to $\mathrm{1}-\exp (-\mathit{\beta }{\mathit{\omega }}_p^q)$, vanishes for infinite temperature.

4 A pseudo-classical version of the Redfield equation

A semi-classical version of the master equation can be extremely useful, allowing one to make use of models derived from the framework of classical mechanics to calculate the spectral density functions. In order to obtain such a theory associated to Eqs. (39) and (42), additional adjustments are necessary. Indeed, because the B^−q(τ) operators do not commute, the correlation functions of the type 〈B^−q(τ)B^q〉 do not obey the general symmetry rules of classical correlation functions. However, symmetrized correlation functions do commute, and so these symmetrized (quantum mechanical) correlation functions should be introduced in order to obtain a semi-classical theory (which we call pseudo-classical to distinguish it from the theory where the effect of the bath is taken into account only through random functions). A general definition of the classical correlation function of two dynamical variables A and B is as follows (Evans and Moriss, 2008):

$$\begin{array}{}\text{(48)}& C_{AB}\left(t\right)=\langle A\left(t\right)B^{\ast }\rangle ,\end{array}$$

where the brackets indicate classical ensemble average. In the case of stationary processes, the following properties of a correlation function can be deduced. Its complex conjugate $C_{AB}^{\ast }\left(t\right)$ is, therefore, as follows (Evans and Moriss, 2008):

$$\begin{array}{}\text{(49)}& \begin{array}{rl}{C}_{AB}^{\ast }\left(t\right)& =\langle A\left(t\right)B^{\ast }{\rangle }^{\ast }=\langle A^{\ast }\left(t\right)B\rangle =\langle A^{\ast }B(-t)\rangle \\ & =C_{BA}(-t).\end{array}\end{array}$$

For an autocorrelation function of A and C_AA(t), one has the following:

$$\begin{array}{}\text{(50)}& C_{AA}(t{)}^{\ast }=C_{AA}(-t).\end{array}$$

Equation (50) shows that, in the general case where the autocorrelation function is complex, $C_{AA}\left(t\right)=C_{AA}^r\left(t\right)+i{C}_{AA}^i\left(t\right)$, with $C_{AA}^r\left(t\right)=\mathit{Re}\left(C_{AA}\right(t\left)\right)$ and $C_{AA}^i\left(t\right)=\mathit{Im}\left(C_{AA}\right(t\left)\right)$ are even and odd functions of time, since $C_{AA}^{\ast }\left(t\right)=C_{AA}^r\left(t\right)-i{C}_{AA}^i\left(t\right)=C_{AA}^r(-t)+i{C}_{AA}^i(-t)$. This implies that the associated spectral density, $J\left(\mathit{\omega }\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{C}_{AA}\left(t\right)e^{-i\mathit{\omega }t}\text{d}t$ is real and $J\left(\mathit{\omega }\right)=J_{AA}^e\left(\mathit{\omega }\right)+J_{AA}^o\left(\mathit{\omega }\right)$, $J_{AA}^e\left(\mathit{\omega }\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{C}_{AA}^r\left(t\right)e^{-i\mathit{\omega }t}\text{d}t$, and $J_{AA}^o\left(\mathit{\omega }\right)=i{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{C}_{AA}^i\left(t\right)e^{-i\mathit{\omega }t}\text{d}t$ are real and, respectively, even and odd functions.

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 follows:

$$\begin{array}{}\text{(51)}& C_{AB}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\langle A^{†}B\left(t\right)\rangle +\langle B^{†}A(-t)\rangle \right\},\end{array}$$

which is real when C_AB(t) is stationary. Note also that the bath operator correlation functions have the following property:

$$\begin{array}{}\text{(52)}& \begin{array}{rl}\langle B^q{B}^{q^{\prime }}\left(\mathit{\tau }\right){\rangle }^{\ast }& =\text{tr}\mathit{\{}{B}^q{e}^{i{H}_{\text{B}}\mathit{\tau }}{B}^{q^{\prime }}{e}^{-i{H}_{\text{B}}\mathit{\tau }}{\mathit{\rho }}^e{\mathit{\}}}^{\ast }\\ & =\text{tr}\mathit{\left\{}{\mathit{\rho }}^e{e}^{i{H}_{\text{B}}\mathit{\tau }}{B}^{q^{\prime }†}{e}^{-i{H}_{\text{B}}\mathit{\tau }}{B}^{q†}\mathit{\right\}}\\ & =\langle B^{q^{\prime }†}\left(\mathit{\tau }\right)B^{q†}\rangle =\langle B^{q^{\prime }†}{B}^{q†}(-\mathit{\tau })\rangle ,\end{array}\end{array}$$

where correlation functions are assumed stationary, so that, in the following:

$$\begin{array}{}\text{(53)}& {\displaystyle }\langle B^q{B}^{-q}\left(\mathit{\tau }\right){\rangle }^{\ast }=\langle B^q\left(\mathit{\tau }\right)B^{-q}\rangle =\langle B^q{B}^{-q}(-\mathit{\tau })\rangle \text{(54)}& {\displaystyle }\langle B^{-q}{B}^q\left(\mathit{\tau }\right){\rangle }^{\ast }=\langle B^{-q}\left(\mathit{\tau }\right)B^q\rangle =\langle B^{-q}{B}^q(-\mathit{\tau })\rangle .\end{array}$$

Using the definitions $C_{\text{L}}^{-q,q}\left(\mathit{\tau }\right)=\langle B^{-q}\left(\mathit{\tau }\right)B^q\rangle$ and $C_{\text{R}}^{q,-q}\left(\mathit{\tau }\right)=\langle B^q{B}^{-q}\left(\mathit{\tau }\right)\rangle$, the relations in Eqs. (53) and (54) show that the average $C^q\left(\mathit{\tau }\right)=\frac{\mathrm{1}}{\mathrm{2}}\left(C_{\text{L}}^{-q,q}\right(\mathit{\tau })+C_{\text{R}}^{q,-q}(\mathit{\tau }\left)\right)$ obey the classical correlation function property $C^{q*}\left(\mathit{\tau }\right)=C^q(-\mathit{\tau })$.

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(τ), as follows:

$$\begin{array}{}\text{(55)}& \begin{array}{rl}{J}^q\left(\mathit{\omega }\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(J_{\text{R}}^{q,-q}\left(\mathit{\omega }\right)+J_{\text{L}}^{-q,q}\left(\mathit{\omega }\right)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(J_{\text{R}}^{q,-q}\left(\mathit{\omega }\right)+e^{-\mathit{\beta }\mathit{\omega }}{J}_{\text{R}}^{q,-q}\left(\mathit{\omega }\right)\right)\\ & ={\displaystyle \frac{\mathrm{1}+e^{-\mathit{\beta }\mathit{\omega }}}{\mathrm{2}}}{J}_{\text{R}}^{q,-q}\left(\mathit{\omega }\right).\end{array}\end{array}$$

Using Eqs. (55) and (45), it is therefore possible to derive an alternative expression of the master equation. This gives the following:

$$\begin{array}{}\text{(56)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _q{S}_p^{-q}{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right){\displaystyle \frac{e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right)\\ & +\sum _q{S}_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right)\\ & +\sum _q{S}^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^{-q}{\displaystyle \frac{e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right)\\ & -\sum _q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{S}_p^{-q}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right),\end{array}\end{array}$$

or the following:

$$\begin{array}{}\text{(57)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& =-\sum _q{S}_p^{-q}{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right){\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right)\\ & +\sum _q{S}_p^{-q}{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right)\\ & +\sum _q{S}_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^{-q}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right)\\ & -\sum _q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q{S}_p^{-q}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}{J}^q\left({\mathit{\omega }}_p^q\right).\end{array}\end{array}$$

Finally, from collecting and rearranging terms, one obtains the following:

$$\begin{array}{}\text{(58)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{J}^q\left({\mathit{\omega }}_p^q\right)\left(\left[S_p^{-q},{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q\right]{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{-\mathit{\beta }{\mathit{\omega }}_p^q}}}\right.\\ & \left.-\left[S_p^{-q},S_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{\mathit{\beta }{\mathit{\omega }}_p^q}}}\right).\end{array}\end{array}$$

In view of clarifying the connection between the derivation of the quantum mechanical master equation to further semi-classical approximations, it is interesting to rewrite Eq. (58) by expanding the temperature function in terms of the parameters $\mathit{\beta }{\mathit{\omega }}_p^q$, as follows:

$$\begin{array}{}\text{(59)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{J}^q\left({\mathit{\omega }}_p^q\right)\\ & \left(\left[S_p^{-q},{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)S_p^q\right]\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathit{\beta }{\mathit{\omega }}_p^q}{\mathrm{4}}}-{\displaystyle \frac{(\mathit{\beta }{\mathit{\omega }}_p^q{)}^{\mathrm{3}}}{\mathrm{48}}}\right)\right.\\ & \left.-\left[S_p^{-q},S_p^q{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\beta }{\mathit{\omega }}_p^q}{\mathrm{4}}}+{\displaystyle \frac{(\mathit{\beta }{\mathit{\omega }}_p^q{)}^{\mathrm{3}}}{\mathrm{48}}}\right)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _q{J}^q\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left[S_p^q,{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right]\right]\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _q{J}^q\left({\mathit{\omega }}_p^q\right)\left({\displaystyle \frac{\mathit{\beta }{\mathit{\omega }}_p^q}{\mathrm{4}}}-{\displaystyle \frac{(\mathit{\beta }{\mathit{\omega }}_p^q{)}^{\mathrm{3}}}{\mathrm{48}}}+\mathrm{\cdots }\right)\\ & \left[S_p^{-q},\left\{S_p^q,{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right)\right\}\right].\end{array}\end{array}$$

Equation (59) contains a double commutator term weighted by the spectral densities $J^q\left({\mathit{\omega }}_p^q\right)$, which are constructed so as to obey the general symmetry properties of classical spectral density functions (see Eq. 48 and below) and are therefore adapted to a semi-classical version of the relaxation master equation. Then, the J^q(ω) are obtained from classical lattice functions of a fluctuating Hamiltonian, whilst the semi-classical master equation obeys detailed balance. The second term of Eq. (59) introduces a lattice-temperature-dependent contribution. However, this term vanishes when the bath operators of the spin–bath coupling Hamiltonian commute, $[B^{-q},B^q]=\mathrm{0}$. According to Eq. (B3), the latter condition also implies that one has the equality of $J_{\text{L}}^{-q,q}\left(\mathit{\omega }\right)=J_{\text{R}}^{q,-q}\left(\mathit{\omega }\right)$, meaning that the lattice temperature is infinite. Stated otherwise, this means that a finite lattice temperature is incompatible with commuting bath operators. However, in general, $[B^{-q},B^q]\ne \mathrm{0}$, so that a detailed balance assumption, or property, which is ensured by the model of a bath in thermal Boltzmann equilibrium, is conveyed to the spin system through the non-commutation of the bath operators B^q. Each term in the series expansion on the right-hand side of Eq. (59) explicitly gives the effect of non-commutation at each order of the parameter $\mathit{\beta }{\mathit{\omega }}_p^q$. The first-order approximation provides the adequate expression in the high temperature limit (see below).

The final relaxation super-operator, which defines the relaxation of the density matrix as ${\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)=\widehat{\mathrm{\Gamma }}{\mathit{\sigma }}^{\ast }\left(t\right)$, may be written as follows:

$$\begin{array}{}\text{(60)}& \widehat{\mathrm{\Gamma }}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _{q,p}{J}^q\left({\mathit{\omega }}_q^p\right){\widehat{\mathcal{D}}}^{\mathit{\beta }}\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},S_p^q\right].\end{array}$$

Here, ${\widehat{\mathcal{D}}}^{\mathit{\beta }}$ is the thermalized double commutator super-operator that composes the super-operator from the two operators A,B in the following way:

$$\begin{array}{}\text{(61)}& {\widehat{\mathcal{D}}}^{\mathit{\beta }}\left(\mathit{\omega }\right)\left[A,B\right]=\left[A,f_{\mathit{\beta }}\left(\mathit{\omega }\right)B•-•B{f}_{\mathit{\beta }}(-\mathit{\omega })\right],\end{array}$$

where $f_{\mathit{\beta }}\left(\mathit{\omega }\right)=\frac{\mathrm{2}}{\mathrm{1}+e^{\mathit{\beta }\mathit{\omega }}}$, and the dot is the place for the operator on which super-operator is applied. It is easy to see that when T→∞, then f_β(ω)→1, and the ${\widehat{\mathcal{D}}}^{\mathit{\beta }}\left(\mathit{\omega }\right)$ becomes a double commutator super-operator, and Eq. (60) becomes a standard sum of the double commutator super-operators. This equation is completely equivalent to Eqs. (42) and (47). Nevertheless, the form of the Lindblad dissipator is still easily recognizable, as one may substitute Eq. (55) into Eq. (42) and obtain the following:

$${\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)=\sum _{p,q}{f}_{\mathit{\beta }}\left({\mathit{\omega }}_p^q\right)J^{q,-q}\left({\mathit{\omega }}_p^q\right)\widehat{\mathcal{L}\mathcal{D}}\left[S_p^{-q},S_p^q\right]{\mathit{\sigma }}_{\text{S}}^{\ast }\left(t\right).$$

Equation (58) partially decouples the statistical and dynamical properties of the heat reservoir. Statistical properties, unrelated to the dynamics, which are functions of the temperature, are contained in the temperature factors, whereas the information about quantum mechanical bath dynamics is contained in the Fourier transform of the symmetrized (here $\frac{\mathrm{1}}{\mathrm{2}}(J_R^{q,-q}+J_L^{-q,q})$) correlation functions. However, the latter still implicitly depend on the temperature through the trace over the bath degrees of freedom.

The case of real correlation functions

The above form of the relaxation master equation (Eqs. 58–60) is suitable for semi-classical approximations of relaxation where classical correlation functions can be used instead of quantum ones that are, in general, impossible to calculate or compute. It is often the case that the classical correlation functions, calculated from classical models of dynamics, such as diffusion, jumps, etc., are real functions of time. The condition of Eq. (50) then implies that the spectral density function J^q(ω) is even $J^q(-\mathit{\omega })=J^q\left(\mathit{\omega }\right)$.

5 Simplifications in the high temperature approximation

When the largest eigenvalue of the operators $S_p^q$ is such that $max\left(\mathit{\beta }{\mathit{\omega }}_p^q\right)\ll \mathrm{1}$, then Eq. (58) takes a simpler form, as follows:

$$\begin{array}{}\text{(62)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _q-J^q\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left[S_p^q,\mathit{\sigma }\left(t\right)\right]\right]\\ & +{\displaystyle \frac{\mathit{\beta }{\mathit{\omega }}_p^q}{\mathrm{2}}}{J}^q\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left\{S_p^q,\mathit{\sigma }\left(t\right)\right\}\right]\end{array}\end{array}$$

where $\mathit{\{}.,.\mathit{\}}$ denotes the anti-commutator.

5.1 The low-order approximation

The assumption that the density operator is always close to the fully disordered state, $\left|\right|\mathit{\sigma }-\frac{\mathrm{1}}{A}\left|\right|\ll \mathrm{1}$, where A is the dimension of the density operator, was made by Redfield (1957). Limiting the expansion to the zeroth order, Eq. (62) becomes the following (Hubbard, 1961):

$$\begin{array}{}\text{(63)}& \begin{array}{rl}{\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _q-J^q\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left[S_p^q,\mathit{\sigma }\left(t\right)\right]\right]\\ & +{\displaystyle \frac{\mathit{\beta }{\mathit{\omega }}_p^q}{A}}{J}^q\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},S_p^q\right].\end{array}\end{array}$$

Moreover, using the property, Eq. (29), and the Taylor expansion of the exponential, it is straightforward to show the following:

$$\begin{array}{}\text{(64)}& e^{-\mathit{\beta }{H}_{\text{S}}}{S}_p^q{e}^{\mathit{\beta }{H}_{\text{S}}}=e^{-\mathit{\beta }{\mathit{\omega }}_p^q}{S}_p^q.\end{array}$$

Therefore, when the density operator is in the thermal equilibrium determined by the Hamiltonian H_S, ${\mathit{\sigma }}^{\mathrm{eq}}={\text{tr}}_{\text{B}}({\exp }^{-\mathit{\beta }{H}_{\text{S}}}{)}^{-\mathrm{1}}{\exp }^{-\mathit{\beta }{H}_{\text{S}}}$, one can show that ℛσ^eq=0, where ℛ is defined by Eq. (58). Then, discarding terms that are second order or higher in $max\left(\mathit{\beta }{\mathit{\omega }}_p^q\right)\ll \mathrm{1}$, in Eq. (62), one obtains the semi-classical formulation of the Redfield equation as follows (Abragam, 1961):

$$\begin{array}{}\text{(65)}& \begin{array}{rl}& {\displaystyle \frac{\text{d}}{\text{d}t}}\left({\mathit{\sigma }}_{\text{S}}^{\ast }\right(t)-{\mathit{\sigma }}^{\mathrm{eq}})=\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{q,p}{J}^q\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left[S_p^q,\mathit{\sigma }\left(t\right)-{\mathit{\sigma }}^{\mathrm{eq}}\right]\right].\end{array}\end{array}$$

The evolution of the expected value of an operator is given by the alternative master equation as follows:

$$\begin{array}{}\text{(66)}& \begin{array}{rl}{\displaystyle \frac{\text{d}}{\text{d}t}}\langle \mathcal{O}\rangle & ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _q-J^q\left({\mathit{\omega }}_p^q\right)〈\left[S_p^q,\left[S_p^{-q},\mathcal{O}\right]\right]〉\\ & +\mathit{\beta }{\mathit{\omega }}_p^q\sum _q{J}^q\left({\mathit{\omega }}_p^q\right)〈\left\{S_p^q,\left[\mathcal{O},S_p^{-q}\right]\right\}〉.\end{array}\end{array}$$

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 the first order in $\mathit{\beta }{\mathit{\omega }}_p^q$. That is, each term in the development is such that, in the following:

$$\begin{array}{}\text{(67)}& \begin{array}{rl}& \text{tr}\left(\left\{S_p^q,\left[\mathcal{O},S_p^{-q}\right]\right\}\mathit{\sigma }\left(t\right)\right)\\ & \ll min\left(\text{tr}\left(\left[S_p^q,\left[S_p^{-q},\mathcal{O}\right]\right]\mathit{\sigma }\left(t\right)\right)\right),\end{array}\end{array}$$

which, at all times, can be discarded. Eq. (67) is, in principle, a less stringent condition and may provide criteria for the quasi- or pseudo-classical approximation, which is 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 obtain an approximate analytical solution.

5.2 The simple case of a two-spin system

5.2.1 Double commutator versus thermal contributions

The differences in the contributions between the first (double commutator) and the second (thermal) series of terms in Eq. (66) are illustrated in Figs. 1 and 2 in the case of a pair of like spins $\frac{\mathrm{1}}{\mathrm{2}}$ subject to the 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, as follows:

$$\begin{array}{}\text{(68)}& J_{\text{DD, ran}}^{\text{cl}}\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{2}{\mathit{\tau }}_{\text{C},\text{ran}}}{\mathrm{1}+{\mathit{\omega }}^{\mathrm{2}}{\mathit{\tau }}_{\text{C, ran}}^{\mathrm{2}}}},\end{array}$$

with correlation times τ_ran=60 ps and τ_C=8 ps for the random field and dipolar interactions, respectively. The dipolar coupling constant $b_{\mathrm{12}}=-({\mathit{\mu }}_{\mathrm{0}}/\mathrm{4}\mathit{\pi }){\mathit{\gamma }}_I^{\mathrm{2}}\mathrm{\hslash }{r}_{\mathrm{12}}^{-\mathrm{3}}$ refers to the dipole–dipole coupling constant, where γ_I is the gyromagnetic ratio, and r₁₂ is the internuclear distance. In these simulations, $b_{\mathrm{12}}=\mathrm{35}\times {\mathrm{10}}^{\mathrm{3}}$ Hz. These values were chosen to give T₁≈2 s and T_S≈20 s.

Figure 1 Expected values of the magnetization spin operator $\mathcal{O}=I_z+S_z$, with the contributions of the double commutator (red curves) and thermal (blue curves) parts of the Redfield relaxation operator to the magnetization from the different operators $S_p^q$ (see Table C1). Panels (a) and (b) correspond, respectively, to the dipolar and random field relaxation of the spins inverted from a Boltzmann equilibrium. Panels (c) and (d) correspond, respectively, to the dipolar and random field relaxation of the spins initially prepared in a singlet state. Simulations were performed using the Scilab software (https://scilab.org, last access: 31 January 2022).

Figure 2 Same as Fig. 1 for the expected value of the singlet operator $\mathcal{O}=\frac{\mathrm{1}}{\mathrm{4}}\mathrm{1}-\mathbf{I}\cdot \mathbf{S}$.

Contributions from both relaxation mechanisms to the expected values in Eq. (66) were computed for the spins prepared either in a singlet state or inverted from thermal equilibrium. The individual terms entering the first and second sums on the right-hand side of Eq. (66) are depicted for the case of the magnetization $\mathcal{O}=I_z+S_z$ (Fig. 1) and singlet state $\mathcal{O}=\frac{\mathrm{1}}{\mathrm{4}}\mathrm{1}-\mathbf{I}\cdot \mathbf{S}$ (Fig. 2) operators. The time evolutions of all the contributions to the rate of change of the expected value of the operator 𝒪(t) are depicted. When the spins are initially prepared in the state $-(I_z+S_z)$, the thermal contribution (blue curve) to the rate of change of the magnetization has no effect, and the only contribution to $\frac{\text{d}〈\mathcal{O}〉\left(t\right)}{\text{d}t}$ comes from the double commutator (red curve). This is the case for both relaxation mechanisms, i.e., dipolar and random field fluctuations (Fig. 1a and b). Moreover, the values chosen for the simulation imply that the dipolar contribution (of the double commutator in this case) to the total relaxation rate $\langle \dot{\mathcal{O}}\rangle \left(t\right)$ is much larger 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 $\langle \dot{\mathcal{O}}\rangle \left(t\right)$ for the dipole–dipole mechanism only (Fig. 1c). In contrast, for terms originating from the random field relaxation, both thermal and dipolar terms are of comparable orders of magnitude (Fig. 1c). These are the terms that cannot be neglected in an approximate solution of Eq. (66). Figure 1c also shows that the dipolar contribution to $\langle \dot{\mathcal{O}}\rangle \left(t\right)$ 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. 1d illustrates the fact that the weight of the thermal contribution decays with time, with the concomitant increase in the double commutator term, which is also due to the progressive depletion of the singlet order.

The situation depicted in Fig. 2 is different and shows the rate of change of the expected value 〈𝒪〉(t), where $\mathcal{O}=\frac{\mathrm{1}}{\mathrm{4}}(\mathbf{I}\cdot \mathbf{S})$, for the same initial state conditions as above. In this case, the dipole–dipole simply does not contribute to $\langle \dot{\mathcal{O}}\rangle \left(t\right)$, as expected from symmetry considerations. This is, of course, the case whatever the initial state (inverted magnetization; see Fig. 2a) or singlet order (Fig. 2c). This illustrates the known fact that singlet state is immune to dipolar relaxation for symmetry reasons.

Alternatively, when the spins are prepared in the $-(I_z+S_z)$ state, both thermal and double commutators contribute, albeit a negligible amount, showing that the spins evolve mostly towards magnetization (compare the scales with Fig. 1b) and that only a negligible part is transferred to a singlet order.

Interestingly, Fig. 2d shows that there is no thermal contribution (blue) to the rate $\langle \dot{\mathcal{O}}\rangle \left(t\right)$, and that, starting from a singlet order, its evolution can be predicted by discarding the thermal terms of Eq. (66) and, therefore, retaining the simple double commutator expression for the relaxation master equation.

5.2.2 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 (1960) procedure, whereby spectral density functions are built from classical ones through the following transformation:

$$\begin{array}{}\text{(69)}& J(\mathit{\omega }{)}_{L,R}\to J^{\text{cl}}(\mathit{\omega })e^{-\frac{\mathit{\beta }\mathit{\omega }}{\mathrm{2}}},\end{array}$$

where J^cl(ω) refers to the classical spectral density function. Equation (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 Bengs and Levitt (2020), J(ω) was assumed to be any kind of spectral density function obtained through classical models, such as diffusion jumps, etc. The distinction between the left and right spectral densities that appear 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 Lindblad'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).

In the following, we derive the evolution of the magnetization of a two-spin system, using the singlet–triplet population basis, and compare the results obtained by both approaches. As above (and in Bengs and Levitt, 2020) the relaxation super-operator $\widehat{\mathrm{\Gamma }}$ is the sum of contributions from mutual dipole–dipole relaxation (${\widehat{\mathrm{\Gamma }}}_{\text{DD}}$) and the interaction with a partially correlated random field (${\widehat{\mathrm{\Gamma }}}_{\text{ran}}$), as follows:

$$\begin{array}{}\text{(70)}& \widehat{\mathrm{\Gamma }}={\widehat{\mathrm{\Gamma }}}_{\text{DD}}+{\widehat{\mathrm{\Gamma }}}_{\text{ran}}.\end{array}$$

The irreducible tensor operator T_λμ representation is better suited to deriving analytical solutions for the problem at hand, where $\widehat{\mathrm{\Gamma }}$ is expressed, according to Eq. (61), as follows:

$$\begin{array}{}\text{(71)}& \begin{array}{rl}& {\widehat{\mathrm{\Gamma }}}_{\text{DD}}={\displaystyle \frac{\mathrm{6}}{\mathrm{5}}}{b}_{\mathrm{12}}^{\mathrm{2}}\sum _{\mathit{\mu }=-\mathrm{2}}^{\mathit{\mu }=\mathrm{2}}{J}_{\text{DD}}^{\text{cl}}\left(\mathit{\mu }{\mathit{\omega }}^{\mathrm{0}}\right){\widehat{\mathcal{D}}}^{\mathit{\beta }}\left(\mathit{\mu }{\mathit{\omega }}^{\mathrm{0}}\right)\left[T_{\mathrm{2}\mathit{\mu }}^{\left(\mathrm{12}\right)†},T_{\mathrm{2}\mathit{\mu }}^{\left(\mathrm{12}\right)}\right],\\ & {\widehat{\mathrm{\Gamma }}}_{\text{ran}}=\sum _{i,j=\mathrm{1}}^{\mathrm{2}}{\mathit{\kappa }}_{ij}{\mathit{\omega }}_{\text{rms}}^{\left(i\right)}{\mathit{\omega }}_{\text{rms}}^{\left(j\right)}\\ & \sum _{\mathit{\mu }=-\mathrm{1}}^{\mathit{\mu }=\mathrm{1}}{J}_{\text{ran}}^{\text{cl}}\left(\mathit{\mu }{\mathit{\omega }}^{\mathrm{0}}\right){\widehat{\mathcal{D}}}^{\mathit{\beta }}\left(\mathit{\mu }{\mathit{\omega }}^{\mathrm{0}}\right)\left[T_{\mathrm{1}\mathit{\mu }}^{\left(i\right)†},T_{\mathrm{1}\mathit{\mu }}^{\left(j\right)}\right].\end{array}\end{array}$$

The T_λμ are eigenoperators of the main Zeeman–Hamiltonian, and their expressions are shown in Appendix C. ${\mathit{\omega }}_{\text{rms}}^{\left(i\right)}$ is the root mean square fluctuation of the random field acting on spin I_i, and in the isotropic case considered here, it is identical for both nuclei, so that ${\mathit{\omega }}_{\text{rms}}^{\left(\mathrm{1}\right)}={\mathit{\omega }}_{\text{rms}}^{\left(\mathrm{2}\right)}$. The coefficient $-\mathrm{1}\le {\mathit{\kappa }}_{\mathrm{12}}\le \mathrm{1}$ describes the degree of correlation of the random field fluctuations on the 1 and 2 nuclei. By definition, ${\mathit{\kappa }}_{\mathrm{11}}={\mathit{\kappa }}_{\mathrm{22}}=\mathrm{1}$. In order to simplify the notations, we will henceforth drop the subscript (κ₁₂→κ).

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_{\mathrm{1}}=R_{\mathrm{1}}^{\text{DD}}+R_{\mathrm{1}}^{\text{ran}}$ are given by, according to Eq. (61), the following:

$$\begin{array}{}\text{(72)}& \begin{array}{rl}{R}_{\mathrm{1}}^{\text{DD}}& =-{\displaystyle \frac{\left(I_z\right|{\widehat{\mathrm{\Gamma }}}_{\text{DD}}\left|I_z\right)}{\left(I_z\right|I_z)}}\\ & ={\displaystyle \frac{\mathrm{3}}{\mathrm{20}}}{b}_{\mathrm{12}}^{\mathrm{2}}{\mathit{\tau }}_{\text{C}}\left(\mathrm{4}{f}_{\mathit{\beta }}\left(\mathrm{2}{\mathit{\omega }}_{\mathrm{0}}\right)+\mathrm{4}{f}_{\mathit{\beta }}(-\mathrm{2}{\mathit{\omega }}_{\mathrm{0}})\right.\\ & \left.+f_{\mathit{\beta }}\left({\mathit{\omega }}_{\mathrm{0}}\right)+f_{\mathit{\beta }}(-{\mathit{\omega }}_{\mathrm{0}})\right)\\ & ={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{b}_{\mathrm{12}}^{\mathrm{2}}{\mathit{\tau }}_{\text{C}},\\ R_{\mathrm{1}}^{\text{ran}}& =-{\displaystyle \frac{\left(I_z\right|{\widehat{\mathrm{\Gamma }}}_{\text{ran}}\left|I_z\right)}{\left(I_z\right|I_z)}}\\ & ={\mathit{\omega }}_{\text{rms}}^{\mathrm{2}}{\mathit{\tau }}_{\text{ran}}\left(f_{\mathit{\beta }}\left({\mathit{\omega }}_{\mathrm{0}}\right)+f_{\mathit{\beta }}(-{\mathit{\omega }}_{\mathrm{0}})\right)\\ & =\mathrm{2}{\mathit{\omega }}_{\text{rms}}^{\mathrm{2}}{\mathit{\tau }}_{\text{ran}},\end{array}\end{array}$$

where $f_{\mathit{\beta }}\left(\mathit{\omega }\right)=\frac{\mathrm{2}}{\mathrm{1}+e^{\mathit{\beta }\mathit{\omega }}}$. Equations (72)–(77) were obtained using the SpinDynamica software (Bengs and Levitt, 2018). It is interesting to note that, in this model, the relaxation rates do not depend on the temperature, which is in contrast to Bengs and Levitt (2020). This is not due to any approximation; rather, it arises from the fact that $f_{\mathit{\beta }}\left(\mathit{\omega }\right)+f_{\mathit{\beta }}(-\mathit{\omega })=\mathrm{2}$, which is explicit in Eq. (72). Similarly, the singlet order relaxation rate is given by the following:

$$\begin{array}{}\text{(73)}& R_{\text{S}}=R_{\text{S}}^{\text{ran}}=-{\displaystyle \frac{\left({\mathbf{I}}_{\mathrm{1}}{\mathbf{I}}_{\mathrm{2}}\right|{\widehat{\mathrm{\Gamma }}}_{\text{ran}}\left|{\mathbf{I}}_{\mathrm{1}}{\mathbf{I}}_{\mathrm{2}}\right)}{\left({\mathbf{I}}_{\mathrm{1}}{\mathbf{I}}_{\mathrm{2}}\right|{\mathbf{I}}_{\mathrm{1}}{\mathbf{I}}_{\mathrm{2}})}}=\mathrm{4}{\mathit{\omega }}_{\text{rms}}^{\mathrm{2}}{\mathit{\tau }}_{\text{ran}}(\mathrm{1}-\mathit{\kappa }),\end{array}$$

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 follows:

$$\begin{array}{}\text{(74)}& \begin{array}{rl}& |S_{\mathrm{0}}\rangle =(|{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\beta }}_{\mathrm{2}}\rangle -|{\mathit{\beta }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}\rangle )/\sqrt{\mathrm{2}},\\ & |T_{+\mathrm{1}}\rangle =|{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}\rangle ,\\ & |T_{\mathrm{0}}\rangle =(|{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\beta }}_{\mathrm{2}}\rangle +|{\mathit{\beta }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}\rangle )/\sqrt{\mathrm{2}},\\ & |T_{-\mathrm{1}}\rangle =|{\mathit{\beta }}_{\mathrm{1}}{\mathit{\beta }}_{\mathrm{2}}\rangle ,\end{array}\end{array}$$

where |α〉 and |β〉 denote the Zeeman spin states of an isolated spin $\mathrm{1}/\mathrm{2}$ (one-half) nucleus with z projection of $+\mathrm{1}/\mathrm{2}$ and $-\mathrm{1}/\mathrm{2}$. In this representation, the population block of relaxation super-operator Eq. (71) is given by Eq. (76).

In Eq. (76), Σ_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 the term $\mathit{\theta }\left(\mathit{\omega }\right)=\text{exp}(-\mathit{\beta }\mathit{\omega }/\mathrm{2})$ by f_β(ω).

In the high temperature limit, both terms become approximately equal, $\mathit{\theta }\left(\mathit{\omega }\right)\approx f_{\mathit{\beta }}\left(\mathit{\omega }\right)\approx \mathrm{1}-\frac{\mathit{\beta }\mathit{\omega }}{\mathrm{2}}$, and the difference between θ(ω) and f_β(ω) becomes significant only in the case of extremely low temperatures or very high frequencies. As could be expected in this limit, the time evolution of the z magnetization is given by the same bi-exponential behavior as in Bengs and Levitt (2020), as follows:

$$\begin{array}{}\text{(77)}& \langle I_z\left(t\right)\rangle /\langle I_z^{\mathrm{eq}}\rangle \approx \mathrm{1}+A_{\mathrm{1}}{e}^{-R_{\mathrm{1}}t}+A_{\text{S}}{e}^{-R_{\text{S}}t},\end{array}$$

with the same coefficients $A_{\mathrm{1}}=\frac{R_{\text{S}}}{\mathrm{2}(R_{\mathrm{1}}-R_{\text{S}})}$ and $A_{\text{S}}=\frac{-\mathrm{2}{R}_{\mathrm{1}}+R_{\text{S}}}{\mathrm{2}(R_{\mathrm{1}}-R_{\text{S}})}$. The fact that A₁ and A_S are the ones found in Bengs and Levitt (2020) is expected because, in this limit, R₁ and R_S do not depend on the temperature factor. In the usual conditions of high but not infinite temperature, it is found that the next nonzero terms in the expansions of A₁(β) and A_S(β) are of the degree β² and, therefore, do not contribute in the regime where ωβ≪1. These expressions can be found in Appendix D.

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, is T≈50 mK. These are unrealistic experimental conditions. Alternatively, the master equation of Eqs. (58)–(61) may well be of use in the context of dynamic nuclear polarization (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 semi-classical theory

In the semi-classical viewpoint (as in 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 time. It is usually understood that this approach does not comply with the Boltzmann equilibrium of the bath. Besides, the stationary state reached by the spin density operator is left undetermined by the master equation so that it must be enforced by the supplementary ad hoc assumption that the spins return to the Boltzmann distribution of the spin populations. A recent analysis by Bengs and Levitt (2020) showed that the usual semi-classical master equation was not able to predict the correct magnetization evolution of a two-spin system prepared in a singlet state.

Thus, the usual semi-classical inhomogeneous master equation provides erroneous predictions in this case. The latter is obtained when the thermal corrections to the double commutator part of the relaxation operator are retained to the first order in the largest eigenvalue $\frac{{\mathit{\omega }}_p^q}{kT}$, and the relaxation operator reduces to a double commutator (low-order case). However, the equilibrium density operator is not a stationary solution in this case, and therefore, a correction term is added to the master equation, leading to the same result as the usual semi-classical master equation (Hubbard, 1961). And so, when the low-order assumption is not verified, as in the case of a spin system prepared in the singlet state, this description becomes inconsistent.

As shown above, the non-commutation of the bath operators has critical consequences, leading to the lattice-temperature-dependent terms in the master equation, and it is only when the bath operators $\left[B^{-q}\right(t),B^q(\mathrm{0}\left)\right]=\mathrm{0}$ that one recovers the double commutator expression, with the additional property J_L(ω)=J_R(ω), so that the Kubo relation imposes an infinite lattice temperature. This illustrates how the finite temperature of the lattice is conveyed to the spins through non-commutation of the bath operators of the coupling Hamiltonian.

The conventional semi-classical approach, where spin–bath interactions are represented by random spin Hamiltonians, has the following two simultaneous consequences: the structure of the relaxation operator is affected in such a way that the master equation takes the form of a double commutator, and since J_L(ω)=J_R(ω), the system cannot evolve to a thermodynamic equilibrium associated with a finite temperature. In this case, the detailed balance property is conserved but only in the special case of infinite lattice temperature. In fact, since the detailed balance is statistical by nature, it is per se compatible with a semi-classical approach. If, on the other hand, a detailed balance is taken into account in the semi-classical theory of NMR relaxation, so that $J(-\mathit{\omega })=e^{-\mathrm{\hslash }\mathit{\omega }/kT}J\left(\mathit{\omega }\right)$ and the general relations Eq. (48) or (50) obeyed by correlation functions are retained, then it is easy to show from the symmetry properties of the spectral density function that the semi-classical Redfield equation (Abragam, 1961) in the following:

$$\begin{array}{}\text{(78)}& {\dot{\mathit{\sigma }}}_{\text{S}}^{\ast }\left(t\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{q,p}J\left({\mathit{\omega }}_p^q\right)\left[S_p^{-q},\left[S_p^q,\mathit{\sigma }\left(t\right)\right]\right],\end{array}$$

is obeyed with this definition of J(ω). However, the expected equilibrium density operator is not a stationary state of Eq. (78) 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 the definition of the correlation function of a phase variable in classical mechanics, $\langle A\left(t\right)B^{\ast }\rangle =\int \text{d}q\text{d}p{\mathit{\rho }}^e{B}^{\ast }{e}^{i\mathbf{L}t}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 (Stillman and Freed, 1980). As stated in several instances in this work, 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. (26) from reducing to the double commutator. Both thermodynamic and quantum mechanical effects are thus entangled in the fully quantum mechanical treatment of relaxation.

Appendix A: Derivation of Eq. (19)

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\text{d}}{\text{d}t}}\mathit{\sigma }\left(t\right)\\ & ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+e^{{\mathcal{L}}_{\text{S}}t}\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{\mathit{\sigma }}^{\ast }\left(t\right)\end{array}\text{(A2)}& {\displaystyle }\begin{array}{rl}& ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)\\ & +\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{e^{{\mathcal{L}}_{\text{S}}t}{\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t\right){\mathcal{L}}_{\mathrm{1}}^{\ast }\left(t^{\prime }\right){\mathit{\rho }}_{\text{B}}^e\right\}{e}^{-t{L}_{\text{S}}}\mathit{\sigma }\left(t\right)\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}& ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{e^{{\mathcal{L}}_{\text{S}}t}{e}^{-({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t}{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t}\right.\\ & \left.e^{-({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t^{\prime }}{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t^{\prime }}{e}^{-t{L}_{\text{S}}}{\mathit{\rho }}_{\text{B}}^e\right\}\mathit{\sigma }\left(t\right)\end{array}\text{(A4)}& {\displaystyle }\begin{array}{rl}& ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{e^{-{\mathcal{L}}_{\text{B}}t}{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t}{e}^{-({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t^{\prime }}\right.\\ & \left.{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t^{\prime }}{e}^{-t{L}_{\text{S}}}{\mathit{\rho }}_{\text{B}}^e\right\}\mathit{\sigma }\left(t\right)\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}& ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t}{e}^{-({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t^{\prime }}\right.\\ & \left.{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})t^{\prime }}{e}^{-t{L}_{\text{S}}}{e}^{-{\mathcal{L}}_{\text{B}}t}{\mathit{\rho }}_{\text{B}}^e\right\}\mathit{\sigma }\left(t\right)\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}& ={\mathcal{L}}_{\text{S}}\mathit{\sigma }\left(t\right)+\underset{\mathrm{0}}{\overset{t}{\int }}\text{d}{t}^{\prime }{\text{tr}}_{\text{B}}\left\{{\mathcal{L}}_{\mathrm{1}}{e}^{-({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})(t^{\prime }-t)}\right.\\ & \left.{\mathcal{L}}_{\mathrm{1}}{e}^{({\mathcal{L}}_{\text{S}}+{\mathcal{L}}_{\text{B}})(t^{\prime }-t)}{\mathit{\rho }}_{\text{B}}^e\right\}\mathit{\sigma }\left(t\right).\end{array}\end{array}$$

​​​​​​​In this derivation, the invariance of the trace to the circular permutations has been used, together with the fact that $[{\mathcal{L}}_{\text{B}},{\mathit{\rho }}_{\text{B}}^e]=\mathrm{0}$, since the bath is in a stationary state.

Appendix B: Evaluation of the terms of Eq. (44)

The correlation functions involved in Eq. (43) are as follows:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\text{R}}^{q,-q}\left({\mathit{\omega }}_p^q\right)\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\langle B^q{B}^{-q}\left(\mathit{\tau }\right)\rangle e^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}\langle f\left|B^q\right|f^{\prime }\rangle \langle f^{\prime }\left|B^{-q}\right|f\rangle \\ & e^{i(f^{\prime }-f)\mathit{\tau }}{e}^{-\mathit{\beta }f}{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^q|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{i(f^{\prime }-f)\mathit{\tau }}{e}^{-\mathit{\beta }f}{e}^{-i{\mathit{\omega }}_p^q\mathit{\tau }}\\ & \approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^q|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i(f^{\prime }-f-{\mathit{\omega }}_p^q)\mathit{\tau }}\text{d}\mathit{\tau }\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^q|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\mathit{\delta }(f^{\prime }-f-{\mathit{\omega }}_p^q)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{L}}\sum _f|\langle f|B^q|f+{\mathit{\omega }}_p^q\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f},\end{array}\end{array}$$

since $B^{-q}=B^{q†}$. Similarly, one has the following:

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_{\text{L}}^{-q,q}\left({\mathit{\omega }}_p^q\right)\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }\langle B^{-q}\left(\mathit{\tau }\right)B^q\rangle e^{-i{\mathit{\omega }}^q\mathit{\tau }}\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}\langle f\left|B^{-q}\right|f^{\prime }\rangle \langle f^{\prime }\left|B^q\right|f\rangle e^{-\mathit{\beta }f}\\ & \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\text{d}\mathit{\tau }{e}^{i(f-f^{\prime })\mathit{\tau }}{e}^{-i{\mathit{\omega }}^q\mathit{\tau }}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^{-q}|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{i(f-f^{\prime }-{\mathit{\omega }}^q)\mathit{\tau }}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^{-q}|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\mathit{\delta }(f-f^{\prime }-{\mathit{\omega }}^q).\end{array}\end{array}$$

Noting that $|\langle f|B^{-q}|f^{\prime }\rangle |=\left|B_{f{f}^{\prime }}^{-q}\right|=\left|B_{f{f}^{\prime }}^{q†}\right|=\left|B_{f^{\prime }f}^{q*}\right|=\left|B_{f^{\prime }f}^q\right|$, and exchanging indices $f\leftrightarrow f^{\prime }$, one obtains, from Eq. (B2) the following:

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^{-q}|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\mathit{\delta }\left(f-f^{\prime }-{\mathit{\omega }}_p^q\right)\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f^{\prime }|B^q|f\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\mathit{\delta }\left(f-f^{\prime }-{\mathit{\omega }}_p^q\right)\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _{f,f^{\prime }}|\langle f|B^q|f^{\prime }\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }{f}^{\prime }}\mathit{\delta }\left(f^{\prime }-f-{\mathit{\omega }}_p^q\right)\\ & ={\displaystyle \frac{\mathrm{1}}{L}}\sum _f|\langle f|B^q|f+{\mathit{\omega }}_p^q\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }\left(f+{\mathit{\omega }}_p^q\right)}\\ & =e^{-\mathit{\beta }{\mathit{\omega }}_p^q}{\displaystyle \frac{\mathrm{1}}{L}}\sum _f|\langle f|B^q|f+{\mathit{\omega }}_p^q\rangle {|}^{\mathrm{2}}{e}^{-\mathit{\beta }f}\\ & =e^{-\mathit{\beta }{\mathit{\omega }}^q}{J}^{\text{R}}\left({\mathit{\omega }}_p^q\right).\end{array}\end{array}$$

Besides, it immediately follows from the definitions of $J_{\text{R}}^q\left(\mathit{\omega }\right)$ and $J_{\text{L}}^q\left(\mathit{\omega }\right)$ so that, in the following:

$$\begin{array}{}\text{(B4)}& J_{\text{R}}^q(-\mathit{\omega })=J_{\text{R}}^{-q}\left(\mathit{\omega }\right).\end{array}$$

Appendix C: Eigenoperators for a homonuclear coupled spin $\mathrm{1}/\mathrm{2}$ pair

In the case of a homonuclear spin pair, the main Hamiltonian is defined as follows:

$$\begin{array}{}\text{(C1)}& H_Z={\mathit{\omega }}^{\mathrm{0}}\left(I_{\mathrm{1}z}+I_{\mathrm{2}z}\right),\end{array}$$

where ${\mathit{\omega }}^{\mathrm{0}}=-\mathit{\gamma }B$, γ is the magnetogyric ratio, and B is the field strength. The eigenoperators for this Hamiltonian are summarized in Table C1. The eigenoperators are denoted by $T_{\mathit{\lambda }\mathit{\mu }}^{\left(ij\right)}$. The superscript (ij) indicates the angular momentum coupling of spins I_i and I_j, resulting in a spherical tensor operator of the total angular momentum λ and z-angular momentum μ.

Table C1 Eigenoperators of Hamiltonian for a homonuclear coupled spin $\mathrm{1}/\mathrm{2}$ pair.

Appendix D: Expansion coefficient for the magnetization evolution solution

The coefficients from the Eq. (77) are in fact temperature dependent, as follows:

$$\begin{array}{}\text{(D1)}& \begin{array}{rl}& A_{\mathrm{1}}\left(\mathit{\beta }\right)=A_{\mathrm{1}}+C_{\mathrm{1}}{\mathit{\beta }}^{\mathrm{2}}+\mathcal{O}\left({\mathit{\beta }}^{\mathrm{4}}\right)\\ & A_{\text{S}}\left(\mathit{\beta }\right)=A_{\text{S}}+C_{\text{S}}{\mathit{\beta }}^{\mathrm{2}}+\mathcal{O}\left({\mathit{\beta }}^{\mathrm{4}}\right),\end{array}\end{array}$$

where the coefficients A₁ and A_S were shown before in the main text, and the coefficients C₁ and C_S are defined as follows:

$$\begin{array}{}\text{(D2)}& \begin{array}{rl}& C_{\mathrm{1}}=\\ & {\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}{\displaystyle \frac{\begin{array}{c}{R}_{\text{S}}\left(\mathrm{600}{R}_{\mathrm{1}}^{\mathrm{3}}-\mathrm{20}{R}_{\mathrm{1}}^{\mathrm{2}}\left(\mathrm{4416}{R}_{\mathrm{1}}^{\text{dd}}+\mathrm{35}{R}_{\text{S}}\right)\right.\\ -\mathrm{5}{R}_{\text{S}}\left(\mathrm{17}\mathrm{280}{\left(R_{\mathrm{1}}^{\text{dd}}\right)}^{\mathrm{2}}+\mathrm{2958}{R}_{\mathrm{1}}^{\text{dd}}{R}_{\text{S}}+\mathrm{5}{R}_{\text{S}}^{\mathrm{2}}\right)\\ \left.+\mathrm{2}{R}_{\mathrm{1}}\left(\mathrm{43}\mathrm{488}{\left(R_{\mathrm{1}}^{\text{dd}}\right)}^{\mathrm{2}}+\mathrm{51}\mathrm{120}{R}_{\mathrm{1}}^{\text{dd}}{R}_{\text{S}}+\mathrm{125}{R}_{\text{S}}^{\mathrm{2}}\right)\right)\end{array}}{\mathrm{5760}{R}_{\mathrm{1}}\left(\mathrm{150}{R}_{\mathrm{1}}-\mathrm{149}{R}_{\mathrm{1}}^{\text{dd}}-\mathrm{25}{R}_{\text{S}}\right){\left(R_{\mathrm{1}}-R_{\text{S}}\right)}^{\mathrm{2}}}}\\ & C_{\text{S}}=-C_{\mathrm{1}}.\end{array}\end{array}$$