Spin relaxation rate for heavy quarks in weakly coupled QCD plasma

We compute the relaxation rate of the spin density of heavy quarks in a perturbative QCD plasma to leading-log order in the coupling constant $g$. The spin relaxation rate $\Gamma_s$ in spin hydrodynamics is shown to be $\Gamma_s\sim g^4\log(1/g)T (T/M)^2$ in the heavy-quark limit $T/M\ll 1$, which is smaller than the relaxation rate of other non-hydrodynamic modes by additional powers of $T/M$. We demonstrate three different methods to evaluate the spin relaxation rate: 1) the Green-Kubo formula in the spin hydrodynamic regime, 2) the spin density correlation function in the strict hydrodynamic limit, and 3) quantum kinetic theory of the spin distribution function in momentum space. We highlight the interesting differences between these methods, while they are ultimately connected to each other by the underlying Ward-Takahashi identity for the non-conserved spin density.


Introduction
The recent experimental observation of the spin polarization of hadrons in the rotating QCD plasma produced in off-central relativistic heavy-ion collisions at RHIC and LHC [1][2][3][4] has motivated much theoretical and phenomenological study of the dynamical evolution of spin polarization in a finite temperature plasma. A large amount of theoretical effort has been invested into the development of a description of spin in kinetic theory [5][6][7][8][9][10][11][12][13][14][15][16][17] and spin in hydrodynamics [18][19][20][21][22][23][24][25][26][27][28][29]. As the spin angular momentum is not conserved and is transferable to orbital angular momentum, the spin density, in general, shows a relaxation behavior toward its local equilibrium value in hydrodynamics [20,28]. The local equilibrium value of spin density is dictated by thermodynamics with rotation, i.e., with conserved total angular momentum [30][31][32][33]. An interesting aspect of this equilibrium state is that the chemical potential corresponding to conserved total angular momentum is equal to the local thermal vorticity of fluids [30,31]. In the strict hydrodynamic regime, the local spin polarization is determined by the local thermal vorticity via a thermodynamic relation. Moreover, spin hydrodynamics in this regime is shown to be equivalent, by a pseudo-gauge transformation [19,[34][35][36], to conventional hydrodynamics with a symmetric energy-momentum tensor with a set of non-dissipative transport coefficients involving fluid vorticity [21,24].
The relaxation rate of the spin density toward its equilibrium value is an important quantity in determining how the spin polarization evolves in time in theoretical simulations of QCD plasma. The QCD plasma produced in heavy-ion collisions goes through several different phases in its lifetime, and one would need to study the dynamics of spin in all -1 -phases to reliably predict the observed value for the spin polarization of hadrons. In this work, we study the spin relaxation rate in the quark-gluon plasma phase, where we assume a high enough temperature to apply a weakly coupled description of QCD plasma, i.e., the finite temperature field theory of perturbative QCD (pQCD) [37][38][39]. Although this is not a realistic assumption for the plasma in heavy-ion collisions where the coupling constant may not be small, the result provides a valuable benchmark in one extreme limit of the theory.
We will focus on the relaxation of spins carried by heavy quarks in the limit M T , where M is the heavy-quark mass, and T is the temperature of the plasma. There are two simplifications in this case. 1) The density of heavy quarks in the plasma is dilute, and we may neglect interactions between heavy quarks. The relaxation of heavy quark spin results from its interactions with other background thermal particles. As we will see, the dominant contribution to the leading-log result comes from scatterings with light hard particles of momentum of order T . 2) The spin relaxation rate in this limit will be shown to be of the order Γ s ∼ g 4 log(1/g)T (T /M ) 2 , where g is the QCD coupling constant. It is smaller than the relaxation rates Γ of other non-hydrodynamic modes by additional powers of T /M 1. This is an important fact that allows us to introduce the spin hydrodynamics in the regime of frequency scale Γ s ω Γ, where the spin density appears as an additional independent quasi-hydrodynamic degree of freedom with a relaxation behavior: an example of the Hydro+ description [40]. The relaxation of spin density in this regime of spin hydrodynamics is governed by a new kinetic coefficient λ s 1 , which appears in the constitutive relations of spin hydrodynamics.
The suppression of spin relaxation rate in the heavy-quark limit can be understood in the non-relativistic limit of heavy-quark dynamics since, in the hydrodynamic limit of a near thermal equilibrium state, the typical heavy-quark velocity is small, v = p/M ∼ T /M 1. In the heavy-quark limit, the QCD Lagrangian involving the heavy quark and gluon is reduced to where ψ is the non-relativistic two-component spinor, D µ ψ = ∂ µ ψ − igA µ ψ with gluon field A µ , B a = abc F bc /2 (a = 1, 2, 3 and abc is the Levi-Civita symbol) is the color magnetic field with the field strength tensor , σ a is the Pauli matrix, L gluon = −(1/2) tr F µν F µν is the gluon Lagrangian, and a term −M ψ † ψ is not shown. The first two leading terms, which are counted as of order T , conserve the non-relativistic SU(2) heavy-quark spin symmetry (see, e.g., ref. [41]), and the leading spin-violating interaction is given by the third term -the Pauli term responsible for the coupling between spin and color magnetic fields. The vital point for our discussion is that the Pauli term is suppressed by the additional power of T /M 1, 2 and this term gives rise to the spin relaxation rate. For our computation of the spin relaxation rate to leading order of the 1/M expansion, it is enough to work with this non-relativistic action, which is simpler than the original relativistic theory of quarks as discussed in ref. [9]. We will henceforth take this non-relativistic theory as the starting point in our computation of the spin relaxation rate in pQCD. Note that the thermal field theory of the gluon field A µ , as well as other light species of quarks, is still relativistic.
There have been a few previous works on the spin relaxation rate in finite temperature pQCD, based on the kinetic-theory framework [9][10][11][12][14][15][16][17] (see also ref. [42] for a recent review of quantum kinetic theory). These works describe the dynamics of the spin of quarks distributed in phase space as a result of scatterings with other thermal particles. Our objective in this work is to compute the kinetic coefficient responsible for the relaxation of spin density, i.e., the heavy-quark rotational viscosity λ s , in the macroscopic description of spin hydrodynamics in the heavy-quark limit M T . Although it is in principle determined by the microscopic collision term obtained in the previous works, the result has not been available in the literature.
We present three different methods in pQCD to compute the spin relaxation rate of heavy quarks in spin hydrodynamics, Γ s = λ s /χ s , where χ s is the spin susceptibility. It is assuring that all three methods give the same result, but the reason why they all should agree is not at all obvious, at least superficially. It is ultimately a consequence of the underlying Ward-Takahashi identity for spin density, together with the separation of scales between Γ s and other relaxation rates Γ, which allows us a subtle controlled limit of frequency ω in the perturbative evaluation of spin density correlation function.
The organization of the paper is as follows. In section 2, we give a brief review of the retarded spin density correlation function predicted by spin hydrodynamics. We discuss the Ward-Takahashi identity for the non-relativistic heavy-quark spin symmetry, which allows us to identify the quantum operator responsible for the spin relaxation in the macroscopic description. Based on this, in section 3, we describe our three different methods of computing the spin relaxation rate Γ s , in the leading-log order of the QCD coupling constant. In section 4, we present the details of computation in these methods and show that they all lead to the same result. Section 5 is devoted to our summary and discussions. In appendix A, we present a derivation of the quantum kinetic equation for heavy quarks based on the Kadanoff-Baym formalism.

Correlation functions in spin hydrodynamics
Our purpose is to consider the dynamics of spin density attached to the non-relativistic fermion ψ. The crucial point here is that the Lagrangian (1.1) enjoys the approximate heavy quark SU(2) symmetry that acts on ψ as ψ → e iθ a σa ψ. Only the Pauli (third) term in eq. (1.1) breaks such heavy quark symmetry, and one can show the following Ward-Takahashi identity as the equation of motion for the spin density operator: ∂ 0 J 0 a + ∇ · J a = Θ a (a = 1, 2, 3), (2.1) where we introduced the spin current J µ a and the source term as Equation (2.1) indicates that spin for the heavy quark is not conserved due to the Pauli term between the spin and color magnetic field. However, it is crucial to note that the source term Θ a has the factor 1/M , so that it is suppressed in the heavy-quark limit. In other words, one can regard eq. (2.1) as an approximate conservation law if M is large enough.
We can coarse-grain operator equation (2.1) to obtain the relationship between the coarse-grained variables for which we shall use the same notations as for the operators. The existence of hydrodynamic description implies that all coarse-grained operators can be expressed as local functionals of the conserved densities such as J 0 a as well as Θ a using constitutive equations. The form of these equations is strongly constrained by the second law of thermodynamics, which requires that the entropy functional S[J 0 a ] = d 3 xs(J 0 a ) is not decreasing in time. Using coarse-grained eq. (2.1), we can express the local entropy production rate as for which we required the local second law of thermodynamics. This means that both the source Θ a and current J i a must be proportional to and its spatial gradients, where we defined internal spin chemical potential µ a , which is a function of spin density J 0 a , as well as external potential for the spin b a , which could be the fluid thermal vorticity, or the external magnetic field times gyromagnetic ratio, or even the torsion. We note that the external spin potential b a is different from the dynamical magnetic field B a . The second law of thermodynamics in eq. (2.3) thus constrains the constitutive equations to leading order in gradients up to semi-positive-definite kinetic coefficients: Substituting constitutive equations (2.5) into the approximate conservation law (2.1), we obtain hydrodynamic equation of motion for spin density J 0 a : Equilibrium is achieved at a value of J 0 a at which the internal and external spin potentials are equal, i.e., µ a = b a . At this point, the entropy is maximized, according to eq. (2.4) and does not increase, according to eq. (2.3). Linearizing equation of motion (2.6) in small deviations around that equilibrium, we find where we introduced a spin diffusion coefficient D s and spin relaxation rate Γ s as using the spin susceptibility χ s defined as usual by Equation (2.7) describes linear response of the spin density to perturbations of the external spin potential, and we can use it to determine the retarded Green's function for spin density [43]: We then find that the retarded spin density correlator has a pole at the imaginary value of which does not vanish in k → 0 limit. Therefore, in hydrodynamic limit k → 0 the spin density shows relaxational behavior with characteristic time τ s = Γ −1 s . If the spin relaxation time τ s is much longer than other microscopic time scales, then the hydrodynamic regime ω Γ can be split into two subregimes. The strict hydrodynamic regime, i.e., ω is much smaller than any relaxation scales, including ω Γ s , and the regime where Γ s ω Γ. The latter is the so-called Hydro+ regime [40], where a small subset of non-hydrodynamic modes relaxes on a scale comparable to the hydrodynamic time scale. Correspondingly, there are two ways to obtain the spin relaxation rate Γ s = λ s /χ s as we shall explain below. By taking a strict hydrodynamic limit ω → 0 at k = 0 of the spin density correlator (2.10) and using eq. (2.8), we obtain one way to evaluate the spin relaxation rate as [28] lim ω→0 Alternatively, we can take the limit ω Γ, while still maintaining ω Γ s : eq. (2.1), Θ a = ∂ 0 J 0 a , the correlator of the spin density J 0 a can be related to the correlator of the source Θ a by a Ward-Takahashi identity, and we obtain (2.14) Using this identity and the fluctuations-dissipation relation we can rewrite eq. (2.13) in terms of the symmetric correlator of Θ a : which gives another useful formula to evaluate heavy-quark rotational viscosity λ s .

Description of the methods
Before we present our detailed computations of the spin relaxation rate in three different methods in section 4, let us give a brief overview of these methods, which summarizes the main ideas, as well as the differences between them. Based on the previous section, let us start with the retarded Green's function of spin density in zero wavenumber limit k → 0: where we use the shorthand notation for the zero wave number Green's function as G(ω) = G(ω, k = 0). The same correlation function computed in finite temperature pQCD must agree with (3.1) in both spin and strict hydrodynamic regimes in ω Γ ∼ g 4 log(1/g)T . The fact that Γ s ∼ Γ T M 2 Γ allows us to consider two different regimes of ω, i.e., the spin hydrodynamic one ω Γ s and the strict hydrodynamic one ω Γ s .
(1) The Green-Kubo formula in spin hydrodynamic regime: In the spin hydrodynamic regime, we have an expansion in powers of Γ s /ω as The first term can also be regarded as the result of taking the controlled limit we discussed in the previous section, i.e., lim Γs ω Γ G . We see that the series is organized in terms of increasing powers of Γ s , which is equivalent to increasing powers of the coupling constant g. Since Γ s → 0 in the perturbative limit g → 0, the condition for the expansion to work, i.e., ω Γ s , is valid for any fixed non-zero value of ω. This means that the diagrammatic perturbation series of the correlation function in field theory should be oneto-one correspondent to the expansion in eq. (3.2), as long as we keep ω fixed and finite. Especially, matching the first term in eq. (3.2) with the first leading diagrams in naive perturbation theory in g, we are able to compute Γ s in simple perturbative computations. We emphasize that this is possible only because the spin density is not a conserved quantity, and we have a finite non-vanishing relaxation rate Γ s that is perturbative in g. This is in -6 -sharp contrast to the correlation functions for conserved quantities, for which the densitydensity correlation function in k = 0 limit simply vanishes identically for all ω due to the Ward-Takahashi identity.
In the diagrammatic computation, it turns out to be easier to compute the Wightman correlation function G where T /ω comes from the limit of Bose-Einstein distribution n B (ω) in the small frequency limit ω T . Therefore, our discussion above implies that the leading diagram for G in the naive perturbation theory should match to 2T χ s Γ s δ ab /ω 2 in Γ s ω T regime. As we discussed in the previous section, one can, in fact, use the Ward-Takahashi identity to relate the correlation functions of spin density J 0 a with those of source operator Θ a , i.e., ∂ 0 J 0 a = Θ a in k = 0 limit. In frequency space, this implies

12
(ω), which we have also checked diagrammatically for our leading order diagrams. The former involves three diagrams, while the latter turns out to be only one diagram in leading order.
We, therefore, use the latter to compute the spin relaxation rate, by matching with the leading order diagram in the naive perturbation theory in g.
Note that this is the Wightman expression of eq. (2.15). As we show in the next section in detail, we are able to determine Γ s in leading-log of coupling constant based on this formula as where m 2 D = g 2 T 2 (2N c + N F )/6 is the Debye mass squared, and C 2 (R) is the Casimir invariant of the color representation R of the heavy quark, which is C 2 (F ) = (N 2 c −1)/(2N c ) for fundamental representation of SU(N c ).
(2) The spin density correlation function in the strict hydrodynamic regime: On the other hand, in the strict hydrodynamic regime of ω Γ s , which is the true hydrodynamic limit of ω → 0, we have an expansion in ω/Γ s as and the same Γ s can be computed from the leading imaginary part in small frequency limit. In this limit, the imaginary part of the retarded correlation function becomes sensitive to infrared singularities arising from the diverging mean-free path of heavy quarks in g → 0 limit, which is called the pinching singularities [44]. The pinching singularity should be regulated by including the imaginary part of selfenergy in the heavy-quark propagators, i.e., the damping rate, which represents the scatterings with thermal background medium that gives rise to a finite mean-free path. The -7 -proper evaluation of correlation functions in this limit further requires a resummation of infinite ladder diagrams [44][45][46], i.e., the vertex corrections, in addition to the damping rate in the propagators: only after this, the Ward-Takahashi identity (or the conservation law in the case of exact global symmetry) is fulfilled [47,48]. One can also understand the pinching singularity from the Boltzmann equation. In the language of kinetic theory, the damping rate and the vertex correction correspond to the loss and the gain terms, respectively, in the collision terms of the Boltzmann equation. The appearance of Γ s in the denominator, which brings in a non-analyticity of the correlation function in coupling g, is a manifestation of pinching singularities of spin density correlation function that is regulated by a finite relaxation rate of spin density.
The fact that only Γ s , not Γ in general, appears in the spin density correlation function after resuming an infinite number of ladder diagrams is also an interesting difference from the case of conserved charges, where the density-density correlation function at k = 0 is simply zero. The case of conserved charges can be understood by replacing Γ s with an infinitesimal that goes to zero. In this limit, the imaginary part linear in ω seems to have a divergent coefficient 1/ , but the condition of expansion ω is never justified. The proper thing to do is to go back to eq. (3.1) with Γ s replaced by , and the density-density correlation function vanishes in → 0 limit. The transport coefficients for conserved charges, e.g., shear viscosity and conductivity, are computed from the current-current correlation functions instead, which shows the similar non-analytic behavior of 1/Γ ∼ 1/[g 4 T log(1/g)] from resummation of infinite ladder diagrams. The Ward-Takahashi identity for spin density should be responsible for how the spin density correlation function has a similar but different non-analytic behavior of 1/Γ s after resummation, which depends only on the spinviolating interactions that are suppressed by powers of T /M . Another way to look at the difference between ω Γ s and ω Γ s in the diagrammatic evaluation of correlation functions is the following. In the expansion in eq. (3.2), which is organized by the naive perturbation theory in g, each term in the series becomes of the same order when ω ∼ Γ s , and one needs to include all terms to properly evaluate the correlation function in ω Γ s limit. Diagrammatically, this necessitates a summation over an infinite number of diagrams. The relevant subset of diagrams in leading order is captured by the diagrams with pinching singularities in the reorganized perturbation theory with damping rate included in the propagators.
(3) The quantum kinetic theory for spin: Our final method of computing the spin relaxation rate of heavy quarks is the quantum kinetic theory of spin density matrix, developed in ref. [9]. The time evolution of a spin distribution function in momentum space, S(p, t), is described by the quantum Boltzmann equation with collision terms, whereΓ S is the quantum collision operator acting on the spin distribution function (not to be confused with Γ s ). The leading-log expression forΓ S is available for a massive quark with its mass satisfying the condition M gT , and we can take M T limit for our -8 -purpose. The result is organized in powers of T /M , (3.8) and the leading term reproduces the momentum diffusion equation with the heavy-quark drag force known in literature [49]. The leading term conserves the total spin density, i.e., pΓ 3 for any S(p), and the spin relaxation in leading order is given by the next termΓ  S . The zero mode represents the equilibrium distribution of spin polarization in momentum space in leading order, which takes a simple Boltzmann distribution. In other words, the zero mode coefficients S 0 corresponds to the spin density in the spin hydrodynamics. How the zero mode relaxes by the spin violating term,Γ (1) S , gives us the spin relaxation rate in hydrodynamics, Γ s . Writing S(p, t) = S 0 (t)e −βEp , and inserting this leading order expression to the quantum kinetic equation (3.8), we obtain, to the leading order in 1/M , (3.9) We used the fact thatΓ with a scalar operatorΓ (1) S , due to rotational invariance of the collision term. We emphasize that the sub-leading corrections in 1/M to the expression, S(p, t) = S 0 (t)e −βEp , exists in general, but are removed upon integration over p, due to the fact that pΓ (0) S [S(p)] = 0 for any S(p). This is the procedure that identifies Γ s correctly to leading order in 1/M in quantum kinetic theory.

Evaluation of the spin relaxation rate in perturbative QCD
In this section, we present our detailed evaluations of the spin relaxation rate in three different ways, all of which lead to the same result given in (3.5). We first compute Γ s via the Green-Kubo formula of the source correlation and the spin density correlation functions in two different regimes, ω Γ s and ω Γ s , in sections 4.1 and 4.2, respectively. In section 4.3, we evaluate the spin relaxation rate in yet another way, based on the quantum kinetic theory of the spin distribution function in momentum space. (ω → 0), evaluated in the naive perturbation theory, and matching that with 2T χ s Γ s δ ab , where χ s is the spin susceptibility [recall eq. (3.4)].
Let us first summarize the basic ingredients for diagrammatic computation. The realtime propagators for the non-relativistic heavy quarks and the relativistic gluons are S 12 (k) = A Derivation of quantum kinetic equation from Kadanoff-Baym formalism In this appendix, we derive the quantum kinetic equation (4.50)-(4.52) for heavy quarks based on the real-time Kadanoff-Baym formalism [59][60][61]. Our starting point is the Schwinger-Dyson equation for heavy-quark Green's function in the real-time formalism, diagrammatically given by where ⌃ ab (a, b = 1, 2) denotes the heavy-quark self-energy in the 12 basis. Note that ra basis are related to 12 basis as After a little cumbersome computation, we can rewrite the Schwinger-Dyson equation depicted in eq. (A.1) as follows: 12 (S (0) ) 1 to this equation, we obtain the left and right Schwinger-Dyson equations (S (0) ) 1 S 12 = ⌃ ra S 12 + ⌃ 12 S ar , (A.4a) We then take the difference of these two equations. Assuming the two-point function S 12 has a weak dependence on the center-of-mass coordinate, we can rely on the derivative expansion. Moreover, recalling the explicit form of the free propagator (S (0) ) 1 = i@ t + 1 2M r 2 , we find that the difference of the Schwinger-Dyson equations becomes where we introduce the Wigner transform of two-point functions and the Moyal product as [XGH: I change the integral variable from x to r.] theory of spin distribution function in momentum space.
4.1 Method 1: the Green-Kubo formula in ! s As we describe in the previous section, the spin relaxation rate can be computed by the source-source correlation function, G S a S b 12 (! ! 0), evaluated in the naive perturbation theory, and matching that with 2T s s ab , where s is the spin susceptibility. We first summarize the basic ingredients for diagrammatic computation. The real-time propagators for the non-relativistic heavy-quarks and the relativistic gluons are where E k = k 2 /2M is the non-relativistic kinetic energy of heavy-quark, ⇢ µ⌫ (k) is the gluon spectral density, and the Fermi-Dirac and Bose-Einstein distribution functions are n F (k 0 ) = 1/(e k 0 e M + 1) and n B (k 0 ) = 1/(e k 0 1), respectively. For simplicity, we use the notation k µ = (k 0 , k) (thus, k 0 corresponds to ! in previous sections), and omit the trivial color structures. Note that the appearance of fugacity e M in n F is due to the fact that the non-relativistic energy is measured from the rest mass energy M . 2 First of all, we note that the source operator S a = g 2M ✏ abc † B b c already contains the small parameter g/M , and it also carries the gluon field. The leading diagram in the naive perturbation theory is the 2-loop sun-set diagram (see figure below). Since the gluon field in S a carries a non-zero frequency-momentum in general (of order gT as we will see), the two propagators of heavy-quark do not carry a same frequency-momentum, and there is no infra-red singularity of pinching poles in the computation. Therefore, we can safely use the free propagator for heavy-quarks in our leading order evaluation.
Moreover, due to the on-shell conditions imposed on the heavy-quark momenta as can be seen in eq. (4.1), the gluon momentum is constrained to be space-like. The free gluon spectral density has support only on light-like momenta. The leading correction to the gluon spectral density ⇢ µ⌫ (k) for space-like momenta comes from the Hard Thermal Loop (HTL) self-energy, which is well-known in literature [? ]. Our computation is based on these ingredients.
The simple scattering picture is also possible by cutting the sun-set diagram in half. As shown in Fig. ??, the diagram simply represents the t-channel scatterings of on-shell heavy-quarks with hard thermal particles of momenta T in the background plasma. It turns 2 We assume that the (relativistic) heavy-quark chemical potential is zero, which leads to a dilute limit theory of spin distribution function in momentum space.
4.1 Method 1: the Green-Kubo formula in ! s As we describe in the previous section, the spin relaxation rate can be computed by the source-source correlation function, G S a S b 12 (! ! 0), evaluated in the naive perturbation theory, and matching that with 2T s s ab , where s is the spin susceptibility. We first summarize the basic ingredients for diagrammatic computation. The real-time propagators for the non-relativistic heavy-quarks and the relativistic gluons are where E k = k 2 /2M is the non-relativistic kinetic energy of heavy-quark, ⇢ µ⌫ (k) is the gluon spectral density, and the Fermi-Dirac and Bose-Einstein distribution functions are n F (k 0 ) = 1/(e k 0 e M + 1) and n B (k 0 ) = 1/(e k 0 1), respectively. For simplicity, we use the notation k µ = (k 0 , k) (thus, k 0 corresponds to ! in previous sections), and omit the trivial color structures. Note that the appearance of fugacity e M in n F is due to the fact that the non-relativistic energy is measured from the rest mass energy M . 2 First of all, we note that the source operator S a = g 2M ✏ abc † B b c already contains the small parameter g/M , and it also carries the gluon field. The leading diagram in the naive perturbation theory is the 2-loop sun-set diagram (see figure below). Since the gluon field in S a carries a non-zero frequency-momentum in general (of order gT as we will see), the two propagators of heavy-quark do not carry a same frequency-momentum, and there is no infra-red singularity of pinching poles in the computation. Therefore, we can safely use the free propagator for heavy-quarks in our leading order evaluation.
Moreover, due to the on-shell conditions imposed on the heavy-quark momenta as can be seen in eq. (4.1), the gluon momentum is constrained to be space-like. The free gluon spectral density has support only on light-like momenta. The leading correction to the gluon spectral density ⇢ µ⌫ (k) for space-like momenta comes from the Hard Thermal Loop (HTL) self-energy, which is well-known in literature [? ]. Our computation is based on these ingredients.
The simple scattering picture is also possible by cutting the sun-set diagram in half. As shown in Fig. ??, the diagram simply represents the t-channel scatterings of on-shell heavy-quarks with hard thermal particles of momenta T in the background plasma. It turns 2 We assume that the (relativistic) heavy-quark chemical potential is zero, which leads to a dilute limit we evaluate the spin relaxation rate in yet another way, based on the quantum kinetic theory of spin distribution function in momentum space.
4.1 Method 1: the Green-Kubo formula in ! s As we describe in the previous section, the spin relaxation rate can be computed by the source-source correlation function, G S a S b 12 (! ! 0), evaluated in the naive perturbation theory, and matching that with 2T s s ab , where s is the spin susceptibility. We first summarize the basic ingredients for diagrammatic computation. The real-time propagators for the non-relativistic heavy-quarks and the relativistic gluons are where E k = k 2 /2M is the non-relativistic kinetic energy of heavy-quark, ⇢ µ⌫ (k) is the gluon spectral density, and the Fermi-Dirac and Bose-Einstein distribution functions are n F (k 0 ) = 1/(e k 0 e M + 1) and n B (k 0 ) = 1/(e k 0 1), respectively. For simplicity, we use the notation k µ = (k 0 , k) (thus, k 0 corresponds to ! in previous sections), and omit the trivial color structures. Note that the appearance of fugacity e M in n F is due to the fact that the non-relativistic energy is measured from the rest mass energy M . 2 First of all, we note that the source operator S a = g 2M ✏ abc † B b c already contains the small parameter g/M , and it also carries the gluon field. The leading diagram in the naive perturbation theory is the 2-loop sun-set diagram (see figure below). Since the gluon field in S a carries a non-zero frequency-momentum in general (of order gT as we will see), the two propagators of heavy-quark do not carry a same frequency-momentum, and there is no infra-red singularity of pinching poles in the computation. Therefore, we can safely use the free propagator for heavy-quarks in our leading order evaluation.
Moreover, due to the on-shell conditions imposed on the heavy-quark momenta as can be seen in eq. (4.1), the gluon momentum is constrained to be space-like. The free gluon spectral density has support only on light-like momenta. The leading correction to the gluon spectral density ⇢ µ⌫ (k) for space-like momenta comes from the Hard Thermal Loop (HTL) self-energy, which is well-known in literature [? ]. Our computation is based on these ingredients.
The simple scattering picture is also possible by cutting the sun-set diagram in half. As shown in Fig. ??, the diagram simply represents the t-channel scatterings of on-shell heavy-quarks with hard thermal particles of momenta T in the background plasma. It turns 2 We assume that the (relativistic) heavy-quark chemical potential is zero, which leads to a dilute limit where E k = k 2 /2M is the non-relativistic kinetic energy of heavy quark, ρ µν (k) is the gluon spectral density, and the Fermi-Dirac and Bose-Einstein distribution functions are n F (k 0 ) = 1/(e βk 0 e βM + 1) and n B (k 0 ) = 1/(e βk 0 − 1), respectively. For simplicity, we use the notation k µ = (k 0 , k) (thus, k 0 corresponds to ω in previous sections), and omit the trivial color structures. Note that the appearance of fugacity e −βM in n F is due to the fact that the non-relativistic energy is measured from the rest mass energy M . 3 First of all, we note that the source operator Θ a = −(g/2M ) abc ψ † B b σ c ψ already contains the small parameter g/M , and it also carries the gluon field. The leading diagram in the naive perturbation theory is the 2-loop sunset diagram [see figure in eq. (4.6) below]. Since the gluon field in Θ a carries a non-zero frequency-momentum in general (of order gT as we will see), the two propagators of heavy quarks do not carry the same frequencymomentum, and there is no infrared singularity of pinching poles in the computation. Therefore, we can safely use the free propagator for heavy quarks in our leading order evaluation.
Moreover, due to the on-shell conditions imposed on the heavy-quark momenta, as can be seen in eqs. (4.1), the gluon momentum is constrained to be spacelike. The free gluon spectral density has support only on light-like momenta. The leading correction to the gluon spectral density ρ µν (k) for spacelike momenta comes from the Hard Thermal Loop (HTL) self-energy, which is well-known in literature [37][38][39]. Our computation is based on these ingredients.
The simple scattering picture is also possible by cutting the sun-set diagram in half. The resulting diagram simply represents the t-channel scattering of an on-shell heavy quark with hard thermal particles of momenta T in the background plasma. It turns out that the leading-log comes from the soft t-channel scatterings, with exchanged gluons carrying momenta in the range m D |k| T . This is similar to other leading-log results for transport coefficients [50,51].
For convenience, we here summarize the HTL results [37] (see, e.g., refs. [38,39] for reviews on the HTL). The HTL spectral density for gluons takes the following form, where the longitudinal and transverse spectral functions are with the corresponding HTL self energies given by Here, m D ∼ gT denotes the Debye mass (more precisely, In the region gT |k| T , and also for the on-shell constraint where v is the heavy-quark velocity of order T /M 1, the spectral functions can be further simplified as [38] which we will use in the following calculation. Based on these ingredients, we now evaluate the spin relaxation rate by our method. Working out the Feynman rule for the source operator involving color magnetic field, we arrive at the following expression for the spin-traced source-source correlation function, reviews on the HTL). The HTL spectral density for gluons takes the following form, where the longitudinal and transverse spectral functions are and ⇢ T (k) = 2 Im with the corresponding HTL self energies given by Here, m D ⇠ gT denotes the Debye mass (more precisely, In the region gT ⌧ |k| ⌧ T , and also for the on-shell constraint where v is the heavy-quark velocity of order p T /M ⌧ 1, the spectral functions can be further simplified as [19] which we will use in the following calculation. Based on these ingredients, we now evaluate the spin relaxation rate by our method. Working out the Feynman rule for the source operator involving color magnetic field, we arrive at the following expression for the spin-traced source-source correlation function,[XGH: It is confusing to use k and k 0 , perhaps to replace k by l in consistence with next subsection.] Using the free propagator for heavy quarks and the HTL spectral density for the gluons, Using the free propagator for heavy quarks and the HTL spectral density for the gluons, we can perform the integral as -11 -where we restricted the |q|-integral in the soft regime discussed above. Note that the q 0 -integral is non-vanishing only in the restricted window, Recalling that the fermion mass M is much larger than any other scales in our problem, we find that the leading-order contribution comes from the q 0 → 0 part of the integrand. Then, using n B (q 0 ) 1/(βq 0 ) and the dilute heavy-quark limit 1 − n F (E p ) 1, and performing both z and q 0 integrals to obtain where we introduced the heavy-quark number density as n(T ) ≡ d 3 p (2π) 3 n F (E p ). Substituting this result into eq. (4.6), we obtain the leading-log result as The spin susceptibility χ s can be computed by introducing a spin chemical potential µ a that couples to the spin in the Hamiltonian as −µ·S(= −µ a S a ) with a spin of heavy quarks S. This induces the modification of distribution functions according to spin direction as n F (E q − µ · S) ≈ n F (E q ) − n F (E q )(µ · S). Then, the induced net spin density along the direction of µ to linear order is 11) where N c is the dimension of the fundamental representation of the color gauge group SU(N c ). As a result, we obtain the spin susceptibility χ s as where we used n F (E) ≈ −βn F (E) in dilute approximation. With this, we finally find the leading-log result of the spin relaxation rate Γ s as where we introduced the Casimir operator for the heavy-quark fundamental representation F as C 2 (F ) = (N 2 c − 1)/(2N c ). By working out color factors with a general representation, it can be shown that the result for a general representation R is given by replacing C 2 (F ) with C 2 (R), leading to the result presented in eq. (3.5).
-12 -4.2 Method 2: Spin density correlation function in ω Γ s regime As explained in section 3, a diagrammatic evaluation of the spin density correlation function in the strict hydrodynamic regime ω Γ s requires a summation of an infinite number of ladder diagrams that are enhanced by infrared pinching singularities. The summation is achieved by solving a Bethe-Salpeter (or Schwinger-Dyson) equation for the vertex corrections given by rungs in the ladder diagrams. In each diagram, the propagators should also include the leading imaginary part of the self-energy, which, via the optical theorem, is proportional to the total cross section of the heavy quark interacting with the background plasma, i.e., the relaxation rate. Thus, the computation proceeds in the same way as the procedure for the perturbative evaluation of transport coefficients [44][45][46][47][48], except for the additional heavy-quark expansion. In the following, we divide our computation into three parts: we first derive the resummed formula for the spin density correlation, next find the Bethe-Salpeter equation for the effective vertex, and finally solve it in the heavy-quark limit.
1) Resummed formula for spin density correlation: We choose to work in the realtime formalism in Schwinger-Keldysh contours [52,53], especially in r/a basis defined by in terms of the doubled fields O 1 and O 2 living on forward and backward contours in time, 1 and 2, respectively. The retarded two point correlation function is related to ra correlation function by G R = iG ra . 4 We would like to compute ra correlation function of spin density, G J 0 a J 0 b ra ≡ J 0 a,r J 0 b,a , in the strict hydrodynamic regime of ω → 0, which is expected to behave as follows [recall eq. (3.6)]: (4.14) The pinchingd singularities affect only the real part of this ra correlation function, that is, the imaginary part of G J 0 a J 0 b R , which contains the spin relaxation rate Γ s . We also note that the imaginary part of G J 0 a J 0 b R corresponds to the spectral function of the spin density thanks to the fluctuation-dissipation relation.
Due to the presence of the pinching singularities, we now need to use the following resummed heavy-quark propagators in ra basis in frequency-momentum space: S ra (k) = We will use the following heavy-quark propagators in ra basis in frequency-momentum space: where ⇣ k ⇠ g 2 log(1/g) denotes the leading relaxation rate for heavy-quark of momentum k, and E k = k 2 2M is the non-relativistic energy. Note that the propagators are in fact 2 ⇥ 2 matrices in spin space, and since they are proportional to the identity matrix in non-relativistic limit, we omit it for simplicity. Once again, the Fermi-Dirac distribution at zero chemical potential in non-relativistic limit is n F (k 0 ) = 1 e k 0 e M +1 ⇡ e M e k 0 ⌧ 1. Recalling J 0a = 1 2 † a , it is straightforward to draw and express real-time Feynman diagrams in the ra basis. The summation over all ladder diagrams is achieved by solving the Bethe-Salpeter equation for the effective vertex function for J 0a , which we denote as 1 2 ⌃ a (q; k), where k is the overall momentum carried by J 0a , and q is the momentum of the out-going propagator from the vertex (see the figure below). Once the effective vertex is found by the Bethe-Salpeter equation, the correlation function G ab ra is given by the following We will use the following heavy-quark propagators in ra basis in frequency-momentum space: where ⇣ k ⇠ g 2 log(1/g) denotes the leading relaxation rate for heavy-quark of momentum k, and E k = k 2 2M is the non-relativistic energy. Note that the propagators are in fact 2 ⇥ 2 matrices in spin space, and since they are proportional to the identity matrix in non-relativistic limit, we omit it for simplicity. Once again, the Fermi-Dirac distribution at zero chemical potential in non-relativistic limit is n F (k 0 ) = 1 e k 0 e M +1 ⇡ e M e k 0 ⌧ 1. Recalling J 0a = 1 2 † a , it is straightforward to draw and express real-time Feynman diagrams in the ra basis. The summation over all ladder diagrams is achieved by solving the Bethe-Salpeter equation for the effective vertex function for J 0a , which we denote as 1 2 ⌃ a (q; k), where k is the overall momentum carried by J 0a , and q is the momentum of the out-going propagator from the vertex (see the figure below). Once the effective vertex is found by the Bethe-Salpeter equation, the correlation function G ab ra is given by the following We will use the following heavy-quark propagators in ra basis in frequency-momentum space: where ⇣ k ⇠ g 2 log(1/g) denotes the leading relaxation rate for heavy-quark of momentum k, and E k = k 2 2M is the non-relativistic energy. Note that the propagators are in fact 2 ⇥ 2 matrices in spin space, and since they are proportional to the identity matrix in non-relativistic limit, we omit it for simplicity. Once again, the Fermi-Dirac distribution at zero chemical potential in non-relativistic limit is n F (k 0 ) = 1 e k 0 e M +1 ⇡ e M e k 0 ⌧ 1. Recalling J 0a = 1 2 † a , it is straightforward to draw and express real-time Feynman diagrams in the ra basis. The summation over all ladder diagrams is achieved by solving the Bethe-Salpeter equation for the effective vertex function for J 0a , which we denote as 1 2 ⌃ a (q; k), where k is the overall momentum carried by J 0a , and q is the momentum of the out-going propagator from the vertex (see the figure below). Once the effective vertex is found by the Bethe-Salpeter equation, the correlation function G ab ra is given by the following where ζ k ∼ g 2 log(1/g) denotes the leading relaxation rate for heavy quark of momentum k, and E k = k 2 /(2M ) is the non-relativistic energy. Note that the propagators are in fact 2 × 2 matrices in spin space, and since they are proportional to the identity matrix in nonrelativistic limit, which we omit for simplicity. Once again, the Fermi-Dirac distribution at zero chemical potential in non-relativistic limit is n F (k 0 ) = 1 e βk 0 e βM +1 ≈ e −βM e −βk 0 1. Recalling J 0 a = 1 2 ψ † σ a ψ, it is straightforward to draw and express real-time Feynman diagrams in the ra basis [46]. The summation over all ladder diagrams is achieved by solving the Bethe-Salpeter equation for the effective vertex function for J 0 a , which we denote as 1 2 Σ a (p; k), where k is the overall momentum carried by J 0 a , and p is the momentum of the out-going propagator from the vertex (see the figure below). Once the effective vertex is found by solving the Bethe-Salpeter equation, the correlation function G J 0 a J 0 b ra is given by the following two real-time diagrams where ⇣ k ⇠ g 2 log(1/g) denotes the leading relaxation rate for heavy-quark of momentum k, and E k = k 2 /(2M ) is the non-relativistic energy. Note that the propagators are in fact 2 ⇥ 2 matrices in spin space, and since they are proportional to the identity matrix in nonrelativistic limit, which we omit for simplicity. Once again, the Fermi-Dirac distribution at zero chemical potential in non-relativistic limit is n F (k 0 ) = 1 e k 0 e M +1 ⇡ e M e k 0 ⌧ 1. Recalling J 0 a = 1 2 † a , it is straightforward to draw and express real-time Feynman diagrams in the ra basis. The summation over all ladder diagrams is achieved by solving the Bethe-Salpeter equation for the effective vertex function for J 0 a , which we denote as where d R is the dimension of the color representation of the heavy-quark (N c for fundamental representation). We recall that the pinching pole enhanced contributions arise from the combination of S ra and S ar in the parallel rail of the diagram. Therefore, we select only a part of S rr written as a linear superposition of S ra and S ar in eq. (4.15), that produces a pinching singularity, where we take the small k 0 ! 0 expansion together with k ! 0, and define n 0 F (q 0 ) ⌘ dn F (q 0 )/dq 0 . To proceed the last line, we use the pinching pole approximation valid at leading order, and also define the on-shell vertex function ⌃ a (q) ⌘ ⌃ a (q 0 = E q , q; k = 0).
where ⇣ k ⇠ g 2 log(1/g) denotes the leading relaxation rate for heavy-quark of momentum k, and E k = k 2 /(2M ) is the non-relativistic energy. Note that the propagators are in fact 2 ⇥ 2 matrices in spin space, and since they are proportional to the identity matrix in nonrelativistic limit, which we omit for simplicity. Once again, the Fermi-Dirac distribution at zero chemical potential in non-relativistic limit is n F (k 0 ) = 1 e k 0 e M +1 ⇡ e M e k 0 ⌧ 1. Recalling J 0 a = 1 2 † a , it is straightforward to draw and express real-time Feynman diagrams in the ra basis. The summation over all ladder diagrams is achieved by solving the Bethe-Salpeter equation for the effective vertex function for J 0 a , which we denote as where d R is the dimension of the color representation of the heavy-quark (N c for fundamental representation). We recall that the pinching pole enhanced contributions arise from the combination of S ra and S ar in the parallel rail of the diagram. Therefore, we select only a part of S rr written as a linear superposition of S ra and S ar in eq. (4.15), that produces a pinching singularity, where we take the small k 0 ! 0 expansion together with k ! 0, and define n 0 F (q 0 ) ⌘ dn F (q 0 )/dq 0 . To proceed the last line, we use the pinching pole approximation valid at leading order, and also define the on-shell vertex function ⌃ a (q) ⌘ ⌃ a (q 0 = E q , q; k = 0).
where d R is the dimension of the color representation of the heavy quark (N c for fundamental representation). We recall that the pinching pole enhanced contributions arise from the combination of S ra and S ar in the parallel rail of the diagram. Therefore, we select only a part of S rr written as a linear superposition of S ra and S ar in eq. (4.15), that produces a pinching singularity, where we take the small k 0 → 0 expansion together with k → 0, and define n F (p 0 ) ≡ dn F (p 0 )/dp 0 . To proceed the last line, we use the pinching pole approximation valid at leading order, and also define the on-shell vertex function Σ a (p) ≡ Σ a (p 0 = E p , p; k = 0). The rotational invariance dictates that the on-shell vertex function divided by the heavy-quark relaxation rate should take the form φ a (p) ≡ Σ a (p) ζ p = σ a f 1 (|p|) + p a (p · σ)f 2 (|p|) + p a f 3 (|p|)1 2×2 + abc p b σ c f 4 (|p|), (4.19) -14 -in terms of four possible scalar functions f i (|p|) (i = 1, 2, 3, 4) that depend on |p| only.
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex Σ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential.
Inserting the above expression of Σ a (p)/ζ p in eq. (4.17), and replacing p a p b with 1 3 p 2 δ ab after angular integration, we arrive at the expression where F (|p|) ≡ f 1 (|p|) + p 2 3 f 2 (|p|). Our remaining task is to find the function F (|p|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex ⌃ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential.
Inserting the above expression of ⌃ a (p)/⇣ p in eq. (4.17), and replacing p a p b with 1 3 p 2 ab after angular integration, we arrive at the expression where F (|q|) ⌘ f 1 (|q|) + q 2 3 f 2 (|q|). Our remaining task is to find the function F (|q|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.

The equation is drawn diagrammatically as
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex ⌃ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential.
Inserting the above expression of ⌃ a (p)/⇣ p in eq. (4.17), and replacing p a p b with 1 3 p 2 ab after angular integration, we arrive at the expression where F (|q|) ⌘ f 1 (|q|) + q 2 3 f 2 (|q|). Our remaining task is to find the function F (|q|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.

The equation is drawn diagrammatically as
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex ⌃ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential.
Inserting the above expression of ⌃ a (p)/⇣ p in eq. (4.17), and replacing p a p b with 1 3 p 2 ab after angular integration, we arrive at the expression where F (|q|) ⌘ f 1 (|q|) + q 2 3 f 2 (|q|). Our remaining task is to find the function F (|q|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.
Deriving Bethe-Salpeter equation: The Bethe-Salpeter equation for ⌃ a (p) is an integral equation, ⌃ a (p) = a + ⌃ a (p), with the source a and the integral kernel ⌃ a (p).  Again, focusing only on the contributions with pinching singularities, we identify the following three diagrams contributing to the integral kernel ∆Σ a (p): (4.19) in terms of four possible scalar functions f i (|p|) (i = 1, 2, 3, 4) that depend on |q| only.

The equation is drawn diagrammatically as
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex ⌃ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential. Inserting the above expression of ⌃ a (p)/⇣ p in eq. (4.17), and replacing p a p b with 1 3 p 2 ab after angular integration, we arrive at the expression where F (|q|) ⌘ f 1 (|q|) + q 2 3 f 2 (|q|). Our remaining task is to find the function F (|q|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.
Deriving Bethe-Salpeter equation: The Bethe-Salpeter equation for ⌃ a (p) is an integral equation, ⌃ a (p) = a + ⌃ a (p), with the source a and the integral kernel ⌃ a (p).  (4.19) in terms of four possible scalar functions f i (|p|) (i = 1, 2, 3, 4) that depend on |q| only.

The equation is drawn diagrammatically as
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex ⌃ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential. Inserting the above expression of ⌃ a (p)/⇣ p in eq. (4.17), and replacing p a p b with 1 3 p 2 ab after angular integration, we arrive at the expression where F (|q|) ⌘ f 1 (|q|) + q 2 3 f 2 (|q|). Our remaining task is to find the function F (|q|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.
Deriving Bethe-Salpeter equation: The Bethe-Salpeter equation for ⌃ a (p) is an integral equation, ⌃ a (p) = a + ⌃ a (p), with the source a and the integral kernel ⌃ a (p).

The equation is drawn diagrammatically as
Again, focusing only on the contributions with pinching singularities, we identify the following three diagrams contributing to the integral kernel ⌃ a (p):  (4.19) in terms of four possible scalar functions f i (|p|) (i = 1, 2, 3, 4) that depend on |q| only.
Upon momentum integration in eq. (4.17), the last two terms do not contribute to G J 0 a J 0 b ra due to angular integration, but there is another reason to expect that f 3,4 in fact vanish. Under parity transformation, the spin vertex ⌃ a remains the same, while the momentum p flips its sign. Therefore, f 3 and f 4 should be absent in a parity-even plasma without, e.g., axial chemical potential. Inserting the above expression of ⌃ a (p)/⇣ p in eq. (4.17), and replacing p a p b with 1 3 p 2 ab after angular integration, we arrive at the expression where F (|q|) ⌘ f 1 (|q|) + q 2 3 f 2 (|q|). Our remaining task is to find the function F (|q|) by solving the Bethe-Salpeter equation for the summation of the ladder diagrams.
Deriving Bethe-Salpeter equation: The Bethe-Salpeter equation for ⌃ a (p) is an integral equation, ⌃ a (p) = a + ⌃ a (p), with the source a and the integral kernel ⌃ a (p).

The equation is drawn diagrammatically as
Again, focusing only on the contributions with pinching singularities, we identify the following three diagrams contributing to the integral kernel ⌃ a (p): where we denote the interaction vertex of heavy quark with gluon field A α as H α (q). Explicitly, these diagrams give us the expression for the integral kernel as where D αβ (q) = A α (q)A β (−q) is the gluon propagator in thermal equilibrium, and C 2 (R) is the second-order Casimir of the heavy-quark representation. In this expression, the interaction vertex H α and the gluon propagator D αβ do not include color generators. Note that the interaction vertex H α (q) consists of three parts, two of which arise from the first two terms in the Lagrangian (1.1), respecting the SU(2) spin symmetry. They are proportional to the identity matrix 1 2×2 . The first term in eq. (1.1) gives Coulomb interaction with A 0 field, and the second term gives the magnetic interaction proportional to the current, which is further suppressed by the small velocity of heavy quarks, |p|/M = v ∼ T /M 1. We thus only need the leading spin conserving interaction in our computation and neglect the second term in the following. On the other hand, we need to keep the third term in eq. (1.1), which is responsible for spin relaxation in leading order. The interaction vertex is therefore given by where the temporal component of the Levi-Civita symbol is zero: 0bc = 0.
We again collect only the pinching pole contributions that arise from pairing S ra and S ar for ∆Σ a (p). Relying on this pinching pole approximation and recalling D αβ rr (q) = 1 2 + n B (q 0 ) (D αβ ra (q) − D αβ ar (q)), we arrive at × D αβ ra (q) − D αβ ar (q) S ar (p + q)S ra (p + q) where D αβ ra (q) − D αβ ar (q) = ρ αβ (q) is the gluon spectral density, for which we use the HTL result summarized in eqs. (4.2)-(4.5) for our leading-log computation. We also use the on-shell condition p 0 = E p , that is imposed in the evaluation of G J 0 a J 0 b ra in eq. (4.17). In terms of the function φ a (p) ≡ Σ a (p)/ζ p , which can be shown to correspond to a linearized distribution function in kinetic theory, the integral equation depicted in eq. (4.21) then takes the following form The important observation, that is common in the similar types of integral equations appearing in the computation of transport coefficients, is that the relaxation rate ζ p from the imaginary part of one-loop self-energy diagram is computed by a very similar expression to the integral kernel, where the factor 1/2 in front comes from the average over spin states since the relaxation rate is independent of the spin states due to rotational invariance. This allows us to combine the two terms in the right-hand side of the integral equation (4.26), and we finally arrive at the Bethe-Salpeter equation for φ a as In the following, we will evaluate this integral equation in leading-log approximation by expanding in powers of small q ∼ gT to quadratic order to derive a second-order differential equation for φ a (p) in p space. This can be considered as the linearized collision term in quantum kinetic theory, perturbed by an external source for the spin density given by the left-hand side of the integral equation. Before moving to the evaluation of Γ s , it is instructive to consider a limit where the SU(2) spin symmetry is respected, by assuming that H α (q) contains only the leading term in 1/M , i.e., the identity matrix in spin space. In this case, the integral equation becomes To see what happens in this limit, we integrate both sides in p-space after multiplying a factor W (p) = n F (E p )[1 − n F (E p )]. By suitable shift of variables and using the properties of the spectral density ρ(q 0 , q) = −ρ(−q 0 , q) = ρ(q 0 , −q), one can show that the right-hand side vanishes after using the identity This parallels the conservation property of total charges in the collision term in kinetic theory. Since W (p) > 0, the left-hand side is not zero after p-integration, which means that the only solution for φ a (p) in this limit is φ a → ∞. This would imply G J 0 a J 0 b ra ∼ χs Γs ωδ ab → ∞, which correctly gives us Γ s → 0 in the limit. From this consideration, it becomes clear that a finite value of Γ s results only from the spin violating interaction in the vertex H α (q) due to the Pauli term. This should be closely related to the (non-)conservative Ward-Takahashi identity for spin density.
3) Solving the Bethe-Salpeter equation in heavy-quark limit: The most laborious part of our computation in this method is the leading-log evaluation of the integral equation (4.28). We first note the presence of a scale hierarchy in the current problem. The leadinglog contributions come from the soft momentum exchange gT |q| T , and the on-shell condition implies q 0 = E p+q − E p ≈ v p · q |q| with v p ≡ p/M , where we use the fact that the typical momentum of heavy quark in eq. (4.17) is p ∼ √ M T from E p ∼ T . This gives us a hierarchy of scales q 0 |q| |p|. In the integral equation, we, therefore, expand each term in powers of |q|/|p| ∼ g, and higher-order terms in q and q 0 bring about terms of higher powers of coupling constant g. To our leading-log computation, it is thus sufficient to expand the terms in the integral equation to second order in q and q 0 [45].
For example, recalling the spin structure of φ a (p) consistent with rotational symmetry, φ a (p) = σ a f 1 (|p|) + p a (p · σ)f 2 (|p|), the expansion of f 1 (|p + q|) would be where θ is the angle between p and q. Similarly, we need for our purpose Furthermore, we can safely neglect n F (E p ) in the following since it is given by n F (E p ) ≈ e −βM e −βEp 1 in the dilute limit of heavy quarks. It is also useful to perform the integration over the angle θ first, by making the energy δ-function into the one for the angle θ, together with restricting the integration range for q 0 as Finally, we use the gluon spectral density in the Coulomb gauge where P ij T (q) = δ ij −q i q j /|q| 2 is the projection operator for spatial components i and j. The HTL resummed longitudinal and transverse spectral densities in the regime q 0 |q| T become the ones given in eq. (4.5), which is sufficient for our leading-log computation.
With all these ingredients used in the integral equation (4.28), performing the leadinglog integral for q and comparing the spin structures in both sides of the equation, we arrive at the two coupled second-order differential equations for f 1 and f 2 as where we define Interestingly, the two equations can be combined to give a single differential equation for that is precisely what we need in order to evaluate the correlation function G -18 -which can be analytically solved in the heavy-quark limit. Before solving this equation in 1/M expansion, it is a useful check to see that the first three terms, that arise from the leading spin-conserving Coulomb interaction, reproduces the conservation property of the collision term in kinetic theory, where the weight factor in our dilute limit becomes W (p) = n F (E p )[1−n F (E p )] ∼ e −βEp = e −p 2 /(2T M ) , up to a constant factor. Therefore, without the last term, which violates the conservation of spin, the solution would diverge F → ∞. This implies that the solution has the following expansion in 1/M as and we are now only interested in the leading solution F (0) . Inserting this expansion into the differential equation (4.38), we obtain the hierarchy of equations It is easy to find the solution F (0) (|p|) = C by inspection, where C is a constant. To determine C, we integrate both sides of the second equation with the weight factor W (p) = e −βEp , that removes the term with F (1) , and we arrive at C = 1/(T 2 Γ). With this leading solution together with eq. (4.20), we finally arrive at the correlation function which should be identified as ω χs Γs δ ab , according to eq. (4.14). Using the spin susceptibility χ s which precisely coincides with eq. (4.13) obtained in the previous subsection.

Method 3: Quantum kinetic theory for spin in heavy-quark limit
Kinetic theory is an intuitive and powerful framework for describing time-dependent dynamics of weakly interacting quasi-particles in phase space (position-momentum space) [54]. It is a statistical description of the system in terms of particle number distribution in phase space. When the two spin states are nearly degenerate in energy, and the quantum correlation time between two spin states is comparable to a macroscopic time scale of interests, one should instead consider the full 2 × 2 density matrix distribution in spin space,ρ(x, p, t), to properly describe time-dependent dynamics of spin polarization of quasi-particles, that is, the quantum kinetic theory. The 2 × 2 density matrix in spin basis can be decomposed into the formρ where f (x, p, t) = tr[ρ(x, p, t)] is identified as the usual particle number distribution function, while S(x, p, t) = tr ρ(x, p, t) σ 2 has the interpretation of the spin distribution function. The time evolution equation for each of them may be written in a similar form as the Boltzmann equation with proper collision terms. Alternatively, the equation written in terms ofρ itself should take a form of the Lindblad equation, in general.
In ref. [9], the quantum Boltzmann equation of f (p, t) and S(p, t) for massive quarks has been constructed to leading-log order of QCD coupling constant, when they are spatially homogeneous, where the collision terms, or the relaxation operators,Γ f andΓ S , are universally of order They take a form of second order differential operator in momentum space p. TheΓ f is nothing but the usual collision term for quark number distribution function in QCD. The collision term for the spin polarization, i.e.Γ S , is a novel object, whose explicit form in spatially homogeneous limit can be found in ref. [9]. We are interested in the heavy-quark limit, i.e., M T , ofΓ S , that determines the relaxation dynamics of heavy-quark spin in momentum space, which ultimately gives us the spin relaxation rate Γ s in the spin hydrodynamic regime. Although it is possible to take M T limit of the result in ref. [9], we present in this section the derivation of Γ S directly from our non-relativistic effective theory (see also the derivation based on the Kadanoff-Baym formalism in appendix A). For clarity of the presentation, we will keep only the leading spin-conserving Coulomb interaction and the leading spin-violating Pauliinteraction, as in the previous section. This is sufficient to obtain the correct spin relaxation rate at the leading order. For more details of our method of derivation, we refer to ref. [9].
The full density operator of heavy quark in spatially homogeneous limit is written in momentum basis asρ whereρ(p, t) is the density matrix in spin space in the subspace of a given momentum p.
More explicitly, we can express it asρ(p, t) = s,s |p, s ρ ss (p, t) p, s |, where |p, s is the -20 - D 22 (t)✓( t) = D 21 (t), and also the fact that the self energy is proportional to identity matrix in spin space due to parity symmetry. One can understand the structure of the above equation (4.50) from the viewpoint of the Lindblad equation [16]: The self energy term corresponds to a dissipative loss of probability, while the cross term is a gain term that ensures the total probability conservation. The index i becomes the tchannel gluon momentum q. One can explicitly check that the total trace of density matrix in momentum-spin space is indeed conserved in our evolution equation. It is also possible to show at this stage that the Boltzmann equilibrium density matrix,⇢ eq (p) / 1 2⇥2 e Ep , is a zero mode of the collision operator. The computation of the above cross and self energy in leading log order is very similar to that in the previous section. Since we are interested in the spin polarization of the density matrix, we focus only on the spin part,⇢(p) ⇠ S(p) · , in the following. Using the gluon spectral density (4.2) in Coulomb gauge, a simple algebra leads to @⇢(p, t) @t where we introduced Note that the C L comes from the Coulomb interaction respecting the spin symmetry, while C T does from the Pauli interaction that breaks it. As in the previous sections, we perform the angular integration of cos ✓ between p and q first, using the relation cos ✓ = M |p| q 0 |q| + |q| 2|p| . This replaces the q integration as -21 - where · · · A denotes the thermal average of background gluon fields interacting with the heavy quark, and the subscript 1 and 2 refers to as the forward and backward Schwinger-Keldysh contours, respectively. This is because the forward contour 1 describes the evolution of ket state |p, s , while the backward contour 2 is for the conjugate bra state p, s |. Note also that the gluon fields in the evolution operatorÛ 1,2 are the Schwinger-Keldysh fields A 1,2 µ on the contours 1 and 2, respectively, whose correlation functions satisfy the fluctuation-dissipation relations in thermal equilibrium.
ExpandingÛ (t+∆t, t) up to the second-order of the quark-gluon interactions in the interaction picture and performing a thermal average of gluon two-point correlation functions, one obtains the time evolution equation for the density matrix as ∂ρ(p, t) ∂t = g 2 C 2 (R)(Γ cross + Γ self energy ), (4.50) where Γ cross arises from the correlation of first order terms in quark-gluon interaction in contours 1 and 2, while Γ self energy is from the second order terms in each contour separately. The diagrammatic representation is depicted in Figure 1. Explicitly, we find where D αβ 12 (q) and D αβ 21 (q) are the thermal gluon correlation functions on Schwinger-Keldysh contours defined in eq. (4.1), and H α (q) is the heavy quark-gluon interaction vertex given in eq. (4.24). In deriving the self-energy term, we need to use the identity D 11 (t)θ(t) + D 22 (t)θ(−t) = D 21 (t), and also the fact that the self-energy is proportional to identity matrix in spin space due to parity symmetry.
One can understand the structure of the above equation (4.50) from the viewpoint of the Lindblad equation [55]: with the anti-commutator {A, B} + ≡ AB + BA. The self-energy term corresponds to a dissipative loss of probability, while the cross term is a gain term that ensures the total probability conservation. The index i becomes the t-channel gluon momentum q. One can explicitly check that the total trace of the density matrix in momentum-spin space is indeed conserved in our evolution equation. It is also possible to show at this stage that the Boltzmann equilibrium density matrix,ρ eq (p) ∝ 1 2×2 e −βEp , is a zero mode of the collision operator.
The computation of the above Γ cross and Γ self energy in leading-log order is very similar to that in the previous section. Since we are interested in the spin polarization of the density matrix, we focus only on the spin part,ρ(p) ∼ S(p) · σ, in the following. Using the gluon spectral density (4.2) in Coulomb gauge, a simple algebra leads to where we introduced C T = (q · σ) q · S(p − q) + q 2 S(p) · σ n B (q 0 ) + q 2 S(p) · σ. (4.55) Note that the C L comes from the Coulomb interaction respecting the spin symmetry, while C T does from the Pauli interaction that breaks it. As in the previous sections, we perform the angular integration of cos θ between p and q first, using the relation cos θ = M |p| q 0 |q| + |q| 2|p| . This replaces the q integration as In computing the qintegration, we have to use q = q L + q T = |q| cos θp + q T with cos θ given as above, while the q T -integration gives zero after polar angle integration. Similarly, q i q j is replaced by The leading-log arises from the expansion of C L,T up to second order in q ∼ gT p, and performing leading-log integral in |q|, with the range from m D ∼ gT to the upper bound of T , where our expression of HTL gluon spectral density breaks down. S , which conserves the total spin density, turns out to be identical to the leading term inΓ f , which describes the momentum diffusion of heavy quark with the known heavy-quark drag force coefficient η D as [49] This results in the Fokker-Planck equation for the density matrix equivalent to the Langevin equation of heavy-quark Brownian motion, with each component S a (p) playing a role of distribution in momentum space, with κ = 2M T η D from the fluctuation-dissipation relation. Physically, what happens is that the spin attached to the heavy quark simply follows the motion of the heavy quark in momentum space. The same result has also been obtained to this order in ref. [12]. The next leading collision operator,Γ S , which is new and encodes the spin-violating effect, is given by a simple expression from which we can evaluate the spin relaxation rate in the leading-log. The leadingΓ (0) S term determines the relaxation rates of non-hydrodynamic modes of spin distribution in momentum space. In other words, its non-zero eigenvalue corresponds to the relaxation rate of the non-hydrodynamic eigenmode. The only exception is the mode with zero eigenvalue, which corresponds to the quasi-hydrodynamic (or Hydro+) mode of spin density. It is easy to guess what the zero mode should be: it is the equilibrium Boltzmann distribution, S 0 e −βEp , with any constant vector S 0 . Indeed, it is a simple algebra to see thatΓ (0) S [S 0 e −βEp ] = 0, while the next order termΓ (1) S gives the relaxation dynamics of spin that we are interested in.
After integration over p, the spin density in position space is given by and S 0 can be identified as the quasi-hydrodynamic mode of spin density. Therefore, we are led to write down the spin distribution function in the 1/M expansion as S(p, t) = S 0 (t)e −βEp + ∆S(p, t), (4.63) where ∆S(p, t) contains all non-hydrodynamic modes which relax much faster than S 0 (t).
It is defined by requiring S(p, t) and S 0 (t) to satisfy the matching condition (4.62), which -23 -means that p ∆S(p, t) = 0. It is expected that ∆S(p, t) is smaller than the leading term by the hydrodynamic expansion parameter ω/Γ. In the regime of spin hydrodynamics where ω ∼ Γ s ∼ (T /M ) 2 Γ, the derivative expansion becomes equivalent to T /M 1 expansion. The relaxation rate of the quasi-hydrodynamic mode, S 0 (t), is the spin relaxation rate Γ s that we are interested in.
After using the expansion (4.63) in the quantum Boltzmann equation (4.47), and integrating over p, we arrive at where we used the matching condition (4.62), as well as the conservation property pΓ (0) S [S(p)] = 0 for any S(p). We dropped the higher-order terms in the 1/M and small frequency expansion, e.g., by neglecting the term such asΓ

Summary and Outlook
In this paper, we have evaluated the spin relaxation rate Γ s for the heavy quark, based on pQCD to leading-log order of coupling constant g. We have formulated three different methods to evaluate Γ s : 1) the Green-Kubo formula based on the source-source correlator in the spin hydrodynamic regime, 2) the spin density correlator in the strict hydrodynamic regime, and 3) the quantum kinetic equation for the heavy-quark spin distribution. While each method demonstrates a different view on spin dynamics, all of these lead to the same result Γ s ∼ g 4 log(1/g)T (T /M ) 2 [see eqs. (4.13), (4.45), and (4.67)]. Thanks to the additional (T /M ) 2 factor, the heavy quark spin shows a parametrically slow dynamics compared to other non-hydrodynamic modes. This scale hierarchy Γ s Γ, where Γ is the relaxation rate for other non-hydrodynamic modes, guarantees the existence of the spin hydrodynamic regime as a well-defined example of Hydro+ [40].
Several outlooks related to the present paper are in order. While we rely on finitetemperature pQCD to evaluate the spin relaxation rate, our formulation -in particular the heavy quark rotational viscosity, λ s in terms of the source or spin correlation functions -is applicable even in the strong-coupling regime, where we have the AdS/CFT correspondence as another theoretical tool [56][57][58][59], and in particular the methods first developed for holographic spin liquids in ref. [60]. It would be interesting to ask how we can formulate heavy-quark spin relaxation in the strong-coupling regime of QCD-like theories, such as N = 4 super Yang-Mills theory, in a manner similar to the heavy quark drag force [61][62][63]. It will give an important benchmark in another extreme limit of the theory, which would be useful in the phenomenological analysis of spin dynamics within the QGP created in heavy-ion collisions.
Another interesting outlook is to apply our result in astrophysics or condensed matter systems. In fact, while we focus on the heavy-quark spin relaxation in QCD plasma in the present paper, our formulation works as well in other systems where spin is approximately conserved. For example, the relaxation rate of proton spin in QED plasma can be studied by our methods with a minor modification. Also, it is interesting to apply the present formulation to various spin systems in condensed matter physics, which deviate from the Heisenberg model only by a small symmetry-breaking perturbation. In this case, the effective field theory with symmetry breaking terms (see, e.g., refs. [64,65]) may be a useful starting point to describe the relaxation rate of spin in condensed matter physics. We leave all these as future work.
x-dependence and the terms with commutators in the left-hand side since the collision term is captured by the right-hand side of this equation. Moreover, expanding the Moyal product simplifies the right-hand side as {A, B} AB + BA. Then, by evaluating the right-hand side (the collision term) with a usual product, we identify this equation as the quantum kinetic equation for heavy quarks including spin, which agrees with eq. (4.50) in the main text.
Let us evaluate the self-energy appearing in the right-hand side of eq. (A.8). As usual, the leading diagram is given by −iΣ ab (p) = @ t + M p · r x S 12 + [i Re ⌃ ra , A ? B + B ?A, respectively. For our purpose, we can neglect all x-dependence and the terms with commutators in the left-hand side since the collision term is captured by the right-hand side of this equation. Moreover, expanding the Moyal product simplifies the right-hand side as {A, B} ? ' AB + BA. Then, by evaluating the right-hand side (the collision term) with a usual product, we identify this equation as the quantum kinetic equation for heavy quarks including spin, which agrees with eq. (4.50) in the main text. Let us evaluate the self-energy appearing in the right-hand side of eq. (A.8). As usual, the leading diagram is given by where the interaction vertex H ↵ (q) between the heavy quark and gluon is given by eq. (4.24).