Open Quantum Dynamics Theory of Spin Relaxation: Application to $\mu$SR and Low-Field NMR Spectroscopies

An open quantum system refers to a system, which is in turn coupled to an environment that can describe time irreversible dynamics through which the system evolves toward the thermal equilibrium state. We present a quantum mechanically rigorous theory in order to help an analysis of spectra obtained from the advanced nuclear magnetic resonance (NMR) and muon spin rotation, relaxation or resonance ($\mu$SR) techniques. Our approach is based on the numerically"exact"hierarchical equations of motion (HEOM) approach, which allows us to study the reduced system dynamics for non-perturbative and non-Markovian system-bath interactions at finite temperature even under strong time-dependent perturbations. We demonstrate the present theory to analyze $\mu$SR and low-field NMR spectra, as an extension of the Kubo-Toyabe theory focusing on the effects of temperature and anisotropy of a local magnetic field, to help further the development of these experimental means.


Introduction
For the analysis of NMR and ESR spectroscopies, the quantum master equation or the Redfield theory has been developed to describe the effects of the longitudinal and transversal relaxations characterized by the time constants T 1 and T 2 . 1,2) Then, the stochastic theory has been employed to describe the effects of the inhomogeneous dephasing characterized by the time decay constant T † 2 in the fast modulation limit. 3) Owing to the advent of experimental techniques that include NMR and ESR, spin dynamics are now investigated under extreme physical conditions, such as quantum computing, where the quantum nature of an environment plays an essential role. 4,5) Thus, such existing theories are insufficient to investigate the complex motion of a spin system. This is also true for zeroto ultralow-field NMR measurement [6][7][8][9] and muon spin rotation, relaxation or resonance (µSR) spectroscopy, 10) because the excitation energy of a spin is almost zero in such measurements and quantum thermalization processes play an important role even at very low temperatures.
µSR spectroscopy is a magnetic resonance technique that utilizes a short-lived elementary particle, a muon (lifetime: 2.2 ×10 −6 s). The muon is a charged spin particle whose magnetic moment is three times larger than that of a proton. Because of its large magnetic moment and short lifetime, it can be implanted in matter to obtain extremely sensitive local magnetic and electronic probes. µSR spectroscopy measures the muon spin polarization recorded from the decay anisotropy of the emitted positrons, as a function of the arrival and decay times. While the experimental setup is completely different from that of NMR spectroscopy, the information obtained by µSR spectroscopy is analogous to that by a low-field NMR measurement.
In 1966, Kubo and Toyabe developed the spin relaxation theory for NMR in zero or weak external magnetic field comparable to the local field from a stochastic approach. 11) Such a low-field measurement was then realized by µSR spectroscopy, and since then, the Kubo-Toyabe theory has been employed to analyze the long-time behavior of the µSR spectrum to probe a local environment of materials. [12][13][14][15][16][17] Various materials that include itinerant helimagnets, superconductors, proteins and DNA have been studied by µSR spectroscopy. [18][19][20][21][22][23][24][25] While several theories for µSR spectroscopy have been developed, [26][27][28][29][30] the Kubo-Toyabe theory is commonly used for investigations of this kind, because it is handy while describing the experimentally obtained µSR signal reasonably well. This feature arises from the assumption that the three-dimensional local random field surrounding the muon is described by a stochastic noise, Ω α (t) for α = x, y, z, which undergoes the Gaussian-Markovian process determined by the noise correlation function, Ω α (t + t 0 )Ω α (t 0 ) = ∆ 2 e −νt , where ∆ and ν are the amplitude and inverse correlation time of the noise, respectively. This allows us to employ the stochastic Liouville equation (SLE) to describe the spin dynamics of the muon. This equation can be solved analytically in a continued fractional form; the static limit of the spin relaxation function is now called the Kubo-Toyabe function.
Although the Kubo-Toyabe theory is convenient to use, there are many limitations in applying it to the analysis of experimental results. For example, this theory does not account for a temperature effect, because the stochastic theory is phenomenological and cannot describe the thermal equilibrium state at finite temperature. It is also applicable only to an isotropic environment without any external forces. Several improvements have been made in the framework of the stochastic theory, 12-17, 29, 30) but applicability is still limited. This is because the stochastic theory relies on the Markovian assumption, whereas the local noise that we investigate arises from the non-Markovian vibrational motion of interand intra-atomic or molecular modes in a complex material.
To eliminate the above-mentioned limitations, here, we consider a system-bath model to treat the system dynamically and use the numerically "exact" hierarchical equations of motion (HEOM) approach to calculate spectra in a rigorous manner. [31][32][33][34][35][36][37][38] The HEOM are the equations of motion that can describe the dynamics of a system for non-perturbative and non-Markovian system-bath interactions at any temperature. In the high temperature limit, the HEOM results for the Drude  bath spectral distribution agree with those from the stochastic theory: The HEOM can be regarded as a generalization of the SLE. Most importantly, the HEOM have flexibility to take into account the effects of a realistic noise that can be obtained from experimental means or molecular dynamics simulations.
This paper is organized as follows. In Sect. 2, we present a typical model system for the NMR and µSR spectroscopy analyses. The HEOM and their characteristic features are described. Numerical results and discussion are presented in Sect. 3. Section 4 is devoted to our conclusions.

Spin-Boson model in three-dimensional space
We consider a spin system as a probe of a local magnetic environment for µSR and NMR spectroscopies described bŷ whereμ ≡ µ 0 sin θ cos φσ x , sin θ sin φσ y , cos θσ z is the dipole operator with the amplitude µ 0 expressed as a function of the solid angle, andσ α (α = x, y, and z) are Pauli matrices. The frequency ω 0 is the Zeeman energy that arises from the static longitudinal external field, which is set to zero in the case of zero-field µSR and NMR spectroscopies. The function B(t) is the time-dependent external field introduced to describe various experimental schemes, which include spin echo, COrrelation SpectroscopY (COSY), and Nuclear Overhauser Effect SpectroscopY (NOESY) measurements utilizing π and/or π/2 pulses. 6) Using B(t), we can explicitly treat nonthermal vibrational motion that is, for example, evaluated from molecular dynamics simulations.
The spin system is independently coupled to three heat baths in the x, y, and z directions to describe an environment in a three-dimensional space (see Fig. 1). We can regard these baths as arising from a local magnetic field owing to the surrounding atoms or molecules. The total Hamiltonian is then given byĤ andĤ α B andĤ α I are the Hamiltonian of the αth bath and the Hamiltonian representing the interaction between the system and the αth bath, respectively. The system part of the systembath interactions is defined asV α =σ α /2, andb α j ,b α † j , ω α j , and g α j are the annihilation operator, creation operator, frequency, and system-bath coupling constant for the jth mode of the αth bath, respectively. For conventional NMR measurements, we consider the heat bath in the z direction only, because the effects of the noise in the x and y directions can be ignored owing to the large ω 0 .
The αth heat bath can be characterized by the spectral distribution function (SDF), defined by By adjusting the form of the SDF, the properties of the local environment consisting of solid-state materials, solvates, and protein molecules can be modeled. The SDF is estimated from spectroscopic experiments [39][40][41] or simulations. [42][43][44][45][46] If we reduce the bath degrees of freedom to obtain the reduced density matrixρ S (t) = tr B {ρ S+B (t)}, the baths produce the noise on the system defined asΩ α (t) ≡ j g α jx α j (t). Through these noise terms, the bath thermalizes the system through fluctuation and dissipation. For a harmonic bath, the effects of thermal fluctuation are expressed as the symmetrized correlation function defined by 31,36) (6) whereas that of the relaxation function is expressed as the anti-symmetrized correlation defined by where · · · B represents the thermal average of the bath degrees of freedom. The symmetrized and anti-symmetrized correlation functions relate through the quantum fluctuationdissipation theorem. The relationship between the present dynamical theory and the stochastic theory can be illustrated using the classical Langevin equation that can be derived from the system-bath model. 36) In the Langevin approach, the func-tionΩ α (t) corresponds to the Langevin random force whose correlation function is defined by Eq. (6). The damping kernel of the Langevin equation is then expressed as Eq. (7), which relates with Eq. (6) through the classical fluctuationdissipation theorem. The stochastic theory thus corresponds to the Langevin formalism without the damping term, because the theory ignores the effects of dissipation. Such a situation is only true when the bath temperature is extremely high and the damping kernel becomes smaller than the fluctuation term.
Because the HEOM formalism treats both fluctuation and dissipation, it can describe the irreversible dynamics of the system accurately, whereas the stochastic theory describes the dephasing motion only. In principle, the HEOM can be constructed for any profile of SDF for the noise correlation functions expressed in terms of damped oscillators as exp[−ζ a ± iω a ], where ζ a and ω a characterize the relaxation and oscillation of the noise correlation for the αth bath, respectively: The HEOM have been derived for the Drude, [31][32][33][34][35][36][37][38] Brownian, [47][48][49] Lorentz, 50) Ohmic, 51) Drude-Lorentz, 52) and their combinations. 53,54) Alternatively, we can derive the HEOM for an arbitrary SDF using the Chebyshev-quadrature spectral decomposition to study the sub-Ohmic SDF at zero temperature. 55) The HEOM for Brownian spectral distribution, which is important to take into account local modes of an environment, are presented in Appendix.

HEOM approach
The HEOM are the equations of motion that allow us to simulate the irreversible dynamics of the system through the fluctuation and dissipation given by Eqs. (9) and (10) in nonperturbative and non-Markovian manners at finite temperature. [31][32][33][34][35][36][37][38] In this formalism, the effects of higher-order non-Markovian system-bath interactions are mapped into the hierarchical elements of the reduced density matrix. This formalism is valuable because it can be used to treat not only strong system-bath coupling but also quantum coherence (quantum entanglement) between the system and the bath, which is essential for studying a system subject to a time-dependent external force and nonlinear response functions. 36) Various analytical and numerical techniques have been developed for the HEOM approaches that allow us to study a complex system under quantum mechanically extreme conditions. With the above described features, the HEOM, which were developed to bridge between the Markovian and perturbative quantum master equation theory and non-Markovian and non-perturbative but phenomenological SLE theory, exhibit wide applicability. The HEOM approach is ideal for extending the applicability of the Kubo-Toyabe low-field theory to various problems and physical conditions in a rigorous manner, and has been applied to spin relaxation problems. 36,[58][59][60] Here, we investigate the µSR problem using the HEOM approach.
In the case of the three-dimensional spin-Boson model, the HEOM is given by [34][35][36][37] ∂ ∂tρ where iL ≡ iĤ × S / and we introduce the set of hierarchy elements n ≡ {n x ; n y ; n z } with n α ≡ {n α 0 , · · · , n α K α } for α = x, y, and z, and the unit vector along the kth element in the α direction expressed as ±e α k that changes the index of the n α k element as n α k ± 1. Here, n α 0 is the element for γ α , whereas n α k for k ≥ 1 are the elements for the Matsubara frequencies ν α k , respectively, in the α direction. The αth bath-induced relaxation operators are defined asΦ where we have introduced the hyperoperator notationÔ ×f ≡ [Ô,f ] andÔ •f ≡ {Ô,f } for any operatorÔ and operand operatorf . The hierarchy of equations of motion introduced above continues to infinity, which is not easy to solve numerically. To truncate Eq.(12), we introduce the terminator 33,34,36) ∂ ∂tρ which is valid for the integers n α 0 , · · · , n α K α satisfying In the high temperature case, the HEOM reduces to 31,36) ∂ ∂tρ where c α 0 in Eq.(13) is now approximated as c α 0 = η α γ α /β and n reduces to n = {n x 0 ; n y 0 ; n z 0 } with e α ≡ e α 0 . Through numerical integration of the equations, we can calculate µSR spectrum under any physical condition even under a time-dependent external force. In the Markovian limit γ α ω c , the above equation further reduces to the master equation In the case of regular NMR described by finite ω 0 , the HEOM can describe the T 1 and T 2 relaxation processes from the x and y baths, respectively, and the T † 2 relaxation process from the z bath in the fast modulation limit without the rotating wave approximation (RWA). While the above equation is valid only in the high temperature case, the HEOM presented in Eq. (12) is valid at any temperature under non-Markovian conditions.
The SLE can also obtained from Eq. (17) by assuming an extremely high temperature case by ignoring the term iγ αV • α /2 and by rescaling ∆ 2 α = η α γ α /β to obtain 31,36) ∂ ∂tρ AAlthough the numerical cost of solving Eq. (19) is almost the same as that of solving Eq. (17), we can solve the above equation in the same manner as Eq. (17) using the truncation scheme developed for the HEOM formalism.

Results and Discussion
We calculated the free induction decay signal of a spin polarization defined by G z (t) = Tr {ρ(t) ·σ z }. To reduce the computational costs, we constructed the HEOM using the Padé-based expression for c α k and ν α k instead of using the Matsubara-frequency-based expression. [61][62][63] Numerical calculations were carried out to integrate Eq. (17) using the fourth-order low-storage Runge-Kutta (LSRK4) method, 64,65) with a time step of δt = 0.01 × 10 −2 . We considered the factorized initial state with the 100% polarized spin in the +z direction, i.e., G z (0) = 1, to account for the condition of the actual µSR measurement. By numerically integrating the SLE presented in Eq. (19), we also calculated the stochastic results to illustrate the roles of the dissipation and the low temperature correction terms involved in the HEOM. We first considered the case of the isotropic environment described by γ x = γ y = γ z = γ and η x = η y = η z = η.

Temperature effects: Interplay between fluctuation and dissipation
In Fig. 2, we show the temperature dependence of the µSR spectrum for the fixed coupling strength η = 0.1 for two cases of the inverse noise correlation time (a) γ = 1 and (b) γ = 10. We compare the HEOM (solid curves) and SLE (dashed curves) results to study the role of the dissipation term by setting ∆ 2 α = η α γ α /β . In the intermediate modulation case shown in Fig. 2(a), the HEOM and SLE results are all similar, whereas in the fast modulation case shown in Fig. 2(b), they are different in the low temperature cases. In the HEOM formalism, the high temperature condition is written as β γ/2 < 1. This implies that all the cases in Fig.  2(a) and the case β = 0.2 in Fig. 2(b) are in the high temperature regime, where the HEOM reduce to Eq. (17). Because we adjusted the amplitude of the stochastic noise to fit the HEOM results, the difference in the HEOM results arises only from the dissipation term presented as the first term in Eq. (13), which becomes negligible for a small β γ in comparison with the second term. This indicates that, by setting ∆ 2 α = η α γ α /β , we may explain the temperature dependence of the µSR spectrum within the framework of the stochastic theory under such conditions. When the temperature becomes very low, however, the signals calculated from the HEOM decay rapidly in comparison with those from the SLE. This difference is due to the time-irreversible dynamics of the spin described by the interplay of the fluctuation and dissipation, whereas the SLE includes dephasing only described by the fluctuation. Because the contribution of the dissipation term becomes large in the low temperature regime, the HEOM results decay more rapidly.

Non-Markovian effects: Role of quantum thermal noise
Next, we investigate the effect of the noise correlation (non-Markovian effects) in the weak coupling case under the (a) intermediate temperature (β = 1.0) and (b) low temperature (β = 2.0) conditions. Because the condition β γ ≤ 0.5 is maintained, the HEOM and SLE results exhibit similar Gaussian decay profiles for the slower modulation case as predicted from the stochastic theory. Note that, when γ is sufficiently small, both HEOM and SLE results exhibit a 1/3 tail that was predicted by the static limit of the Kubo-Toyabe theory. The distinct feature of the HEOM results is observed in the fast modulation cases (γ=10) in Figs. 3(a) and 3(b): The signals calculated from the HEOM decay more slowly than those calculated from the SLE after exhibiting a fast initial decay in the time period less than 1/β ≈ 1.0 or 0.5.
While the fast decay is due to the population relaxation arising from the dissipation, the slow decay is due to the quantum dephasing arising from the quantum thermal noise. As illustrated in Eqs. (9) and (10), two types of non-Markovian noise are involved in the system dynamics: one is of mechanical origin characterized by the fluctuation (c α 0 e −ν α 0 t ) and dissipation (c α 0 e −ν α 0 t ) with ν α 0 = γ, and the other is of quantum thermal origin characterized by the fluctuation only (c α k e −ν α k t for k ≥ 1) with ν α k = 2πk/β . When γ is much larger than ν 1 , the me- chanical contribution with e −ν α 0 t vanishes after t > 1/ν 1 , and the effects from the quantum thermal noise take place. The quantum thermal fluctuation exhibits a peculiar behavior in comparison with the mechanical fluctuation, because the amplitude of the noise becomes negative for a large a γ [see Fig.  7(b) in Ref. 36]. Thus, the signal obtained from the HEOM decays more slowely than that obtained from the SLE. Although the SLE is also a non-Markovian theory, this quantum thermal dephasing process can be described only from the numerically "exact" HEOM approach.

Non-perturbative system-bath interactions
We study the non-perturbative effects of the system-bath coupling by changing η. As illustrated in Fig. 4(a), the differences between the HEOM and SLE results increase with coupling strength even in the intermediate modulation case. This is because, while the amplitudes of the fluctuation and dissipation are both proportional to the coupling strength, the relaxation arising from the dissipation plays a greater role than the dephasing arising from the fluctuation, owing to the timeirreversible nature of the relaxation. As depicted in Fig. 4(b), such differences become prominent in the faster modulation case, as in the cases described in Sects. 3.1 and 3.2. The time period of the initial decay decreases with increasing coupling strength, because the noise with c k e −ν α k t for a larger k can interact with the system several times in this non-perturbative regime.

Anisotropic effects of environment
The HEOM formalism is ideal for studying a spin system under realistic conditions, because it allows the treatment of various anisotropic environments with any profile of noise correlation functions characterized by for any combination of α, α = x, y, z. This is because the HEOM formalism is based on the equations of motion approach. Below, we investigate the anisotropy effects of noise amplitudes and noise correlations. For this purpose, we consider the extremely high temperature case (β = 2.5 × 10 −3 ) with the weak system-bath coupling η = 1.0 × 10 −3 . Thus, the HEOM and SLE results become almost identical for ∆ 2 α = η α γ α /β .

Anisotropic noise amplitudes
We first study the effects of the anisotropic system-bath coupling strength expressed as where a, b, and c are the anisotropic constants. While the analysis of anisotropic effects was limited in the static case on the basis of the Kubo-Toyabe theory, 16,17) there is no technical limitation from the HEOM approach for such problems, because we are only integrating the HEOM.
In Fig. 5 we present the signals for the (i) x anisotropic (a = 1, b = 2, and c = 1), (ii) z anisotropic (a = 1, b = 1, and c = 2), and (iii) isotropic (a = 1, b = 1, and c = 1) cases for various γ values. The other parameters are fixed as β = 2.5 × 10 −3 and η = 1.0 × 10 −3 . In Figs. 5(a) -5(c), the signal decays more rapidly in the x anisotropic case than in the isotropic case, whereas the signal decays more slowly in the z anisotropic case than in the isotropic case. This is because theσ x operator causes the longitudinal (T 1 ) relaxation, whereas theσ z operator causes not a relaxation but a dephasing (T † 2 ) for the z-polarized spin. When γ increases, the spin distribution approaches the equilibrium value owing to the relaxation. In the fast modulation case shown in Fig. 5(d), the signal decays more rapidly than in the case shown in Fig. 5(c), because the effective coupling strength becomes larger for a large γ owing to the factor γ 2 α /(γ 2 α + ω 2 c ), where ω c is the characteristic frequency of the system dynamics. The z anisotropic results become similar to the isotropic case, because when the spin element in the z direction becomes small, the effects of dephasing in the z direction also become minimal.

Anisotropic noise correlation
Finally, we consider the case that some of the threedimensional baths are correlated. Such a model was developed to analyze the noise correlation of different vibrational modes by two-dimensional infrared spectroscopy. 66,67) Here, we consider the (i) x − y correlated [V x−y = (σ x +σ y )/2 and V z =σ z ], and (ii) x − z correlated [V x−z = (σ x +σ z )/2 and V y =σ y ] cases for the bath Hamiltonian, Eq. (4), with (i) α = x − y, and z, and (ii) α = x − z, and y. These results are compared with (iii) the isotropic case α = x, y, and z.
In Fig. 6, µSR spectra are presented for various inverse noise correlation times γ = (a) 0.1, (b) 0.3, (c) 1, and (d) 10. The characteristic feature of the present results is illustrated from the fluctuation term in Eq. (18) expressed as In the x − y correlation case, this term is expressed asV × x +Ĉ x−y +σ ×2 y )/4, whereĈ x−y = (σ × xσ × y +σ × yσ × x ). In the slow modulation case shown in Fig.  6(a), the movements of the spin in the x and the y direc-tions are not correlated, and the contribution fromĈ x−y becomes small. Thus, the longitudinal (T 1 ) relaxation becomes weaker than that in the isotropic case. In the fast modulation case shown in Fig. 6(d), however, the contribution fromĈ x−y becomes similar to that ofσ ×2 x andσ ×2 y , and thus we havê y , which leads to the x − y result becoming similar to the isotropic one. In the x − z correlation case, we haveV × under the slow modulation condition, whereas the contribution from the y direction does not change. Thus the x − z signal shown in Fig. 6(a) decays more slowely than that in the isotropic case, whereas it still decays more rapidly than the x − y correlated signal. Under the fast modulation condition shown in Fig. 6(d), the decay of the signal in the x−z correlated case is slow, because, for the z-polarized spin, theĈ x−z contribution remains small even under the fast modulation condition.
The above results indicate that the µSR spectrum is sensitive to the anisotropic effects of the environment, which should be detected experimentally in accordance with the theoretical analysis.

Conclusions
As illustrated in this paper, the HEOM approach has distinct features for the analysis of µSR and near-zero-field NMR spectra. First, while the stochastic approach can treat the hightemperature Markovian case only, the HEOM approach can treat the realistic non-Markovian noise arising from complex environments, such as nanomaterials, proteins, and a spin lattice in different magnetic ordered phases. This is because the HEOM are constructed on the basis of a fairly complex system-bath Hamiltonian: It is also possible to construct a simulation model on the basis of a molecular dynamics simulation. [42][43][44][45][46] Second, because the HEOM is a dynamical theory, we can easily and clearly identify the roles of the systembath interaction, noise correlation time, and heat-bath temperature. Third, because the HEOM approach is the equations of motion approach, there is no difficulty in taking into account the effects of a time-dependent external field. This allows us to calculate multi-dimensional near zero-field NMR signals for various pulse sequences. In addition, we can include the effects of nonthermal local environmental modes explicitly as the time-dependent external field, whereas the other thermal effects are taken into account using the hierarchical structure. For numerical integration, we can employ a complex quantum system, such as a spin chain 58) or a spin lattice 59, 60) as the main system instead of a simple spin system.
In conclusion, the present formalism provides a powerful means of analyzing µSR and various NMR measurements for the study of environmental effects. It is also possible to use the HEOM theory to investigate other scattering and spectroscopic measurements, which include neutron scattering, electron paramagnetic resonance (EPR), and Mössbauer measurements. 36) All of the possibilities mentioned above can be carried out as future studies upon request.
Here and for all previously defined abbreviations throughout. Thus, we can construct the HEOM by evaluating the time derivative of the reduced density matrices as 48,49) ∂ ∂tρ The hierarchical elements n ≡ {n x ; n y ; n z } are now defined by n α ≡ {n α 0 ,n α 0 , n α 1 , · · · , n α K α }, where n α 0 andn α 0 are the elements for ν α 0 andν α 0 for α = x, y, and z, respectively. The unit vectors in the α direction, which change the indexes of the n α 0 andn α 0 elements as n α 0 ± 1 andn α 0 ± 1, are expressed as ±e α and ±ē α , respectively, whereas the other unit vectors e α k are defined in the same manner as Eq. (12).