Spin-thermal shear coupling in a relativistic fluid

We show that spin polarization of a fermion in a relativistic fluid at local thermodynamic equilibrium can be generated by the symmetric derivative of the four-temperature vector, defined as thermal shear. As a consequence, besides vorticity, acceleration and temperature gradient, also the shear tensor contributes to the polarization of particles in a fluid. This contribution to the spin polarization vector, which is entirely non-dissipative, adds to the well known term proportional to thermal vorticity and may thus have important consequences for the solution of the local polarization puzzles observed in relativistic heavy ion collisions.


I. INTRODUCTION
In a rotating fluid at global thermodynamic equilibrium, particle spin gets polarized along the direction of the angular velocity vector by an amount which is proportional to ω/KT . This phenomenon is the essence of the Barnett effect [1] and it has been known for a long time. In a relativistic fluid at local thermodynamic equilibrium, the covariant form of statistical mechanics dictates that spin polarization is driven by thermal vorticity: where β is the four-temperature vector: u being the four-velocity and T the proper temperature. At first order in thermal vorticity, the formula relating the mean spin vector S µ (p) of a spin 1/2 fermion to thermal vorticity reads [2]: where Σ is a 3D hypersurface, n F is the Fermi-Dirac phase-space distribution function: q being the charge of the particle and µ the corresponding chemical potential. The equation (3) predicts that a particle can get a spin polarization in the presence of gradients of temperature, vorticity and acceleration. The observation of spin polarization in relativistic nuclear collisions [3] confirmed the predictions of the formula (3) for the global polarization (with Σ the hadronization hypersurface), that is integrated over all momenta. In fact, the formula (3) failed to reproduce the measurements as a function of momentum [4][5][6]; particularly, the sign of the longitudinal polarization and the polarization along the angular momentum as a function of the azimuthal angle, which has been investigated in several papers [7][8][9][10][11][12][13][14][15].
In this work, we will show that also the symmetric gradient of β contributes to the spin at local thermodynamic equilibrium at the leading order. This term is non-dissipative as well as non-local for it depends on a specific 3D integration hypersurface and implies a new relativistic effect, namely a coupling between spin and the shear tensor in a relativistic fluid.

Notation
In this paper we adopt the natural units, with = c = K = 1. The Minkowskian metric tensor g is diag(1, −1, −1, −1); for the Levi-Civita symbol we use the convention 0123 = 1. We will use the relativistic notation with repeated indices assumed to be saturated. Operators in Hilbert space will be denoted by a wide upper hat, e.g. H, except the Dirac field operator which is denoted by a Ψ.

II. LOCAL THERMODYNAMIC EQUILIBRIUM AND ITS GRADIENT EXPANSION
For a relativistic quantum fluid which, at some time, achieves Local Thermodynamic Equilibrium (LTE), a powerful approach is the Zubarev's method of the stationary non-equilibrium density operator [16,17]. We refer the reader to the recent paper [18] for a more detailed description. The actual density operator of such a fluid, in the Heisenberg representation, is: where β is the four-temperature vector, ζ the ratio between chemical potential and temperature and Σ eq is some initial 3D hypersurface where LTE is achieved. For relativistic nuclear collisions, this is supposedly the 3D hyperbolic hypersurface where the quark-gluon plasma (QGP) thermalizes (see figure 1). It should be pointed out that the form of the local equilibrium density operator is pseudo-gauge dependent [19,20], with the form in eq. (4) applying to the Belinfante stress-energy tensor only. We note right away that the final result of this work would be the same if we used the canonical stress-energy tensor instead; this is shown in detail in the Appendix A. Henceforth, it will be understood that T is the Belinfante symmetrized stress-energy tensor. The operator (4) can be turned into a more manageable form by means of the Gauss' theorem, taking into account that T and j are conserved currents: where Σ(τ ) is some 3D hypersurface at "present" time τ . In the case of heavy ion collisions the hypersurface Σ(τ ) is usually the joining of the freeze-out hypersurface Σ F O encompassing the QGP space-time region and the two side branches σ ± subsets of the Σ eq , as shown in fig. 1. A peculiarity of the heavy ion collisions is that the hypersurface of "present" local equilibrium is partly time-like, that isn ·n = −1.
In the right hand side of the density operator (5) the first term is the LTE and, as expected for a quasi-ideal fluid such as the QGP, it is the predominant one; the second term is, on the other hand, supposedly a correction and it is responsible for everything that can be called dissipative. Indeed, it can be shown that entropy is generated only if the second term is non-vanishing [16].
We shall focus on the LTE term: The density operator (6), as well as eqs. (4) and (5), is independent of the hypersurface only if the β field satisfies the Killing equation [21]. When calculating the mean values of any local operator O(x), as the thermodynamic fields β and ζ are supposedly slowly varying compared to the correlation lengths between O and the operators T and j (see e.g. ref. [22]), it is a good approximation to expand β and ζ in a Taylor series from x: and similarly for ζ. Thus, we have: For the sake of simplicity, we will omit the gradients of ζ and focus on the gradients of β, which are the most relevant for our purposes. These gradients can be split into a symmetric and an anti-symmetric part giving rise to: where is the thermal vorticity (1). We can recognize in the first term of the above equation the total angular momentum-boost operator J λν x (with a proviso, see ref. [23]) centered in x, while the second term includes the non-conserved operator: coupled to the thermal shear tensor: which vanishes at global thermodynamic equilibrium due to the Killing condition. It is very important to stress that Q x is a tensor in a more limited sense than the angular momentum-boost operator J x . Indeed, since the integrand of Q x in eq. (9) is not divergenceless: its value specifically depends on the hypersurface Σ, unlike J x , and, strictly speaking, should then be denoted as Q x (Σ) (even though we will not use that notation). In quantum language, the operator Q x does not fulfill the transformation rule for a tensor operator under a Lorentz transformation, that is: where Λ is the unitary representation of the Lorentz transformation Λ in the Hilbert space. This implies that all results involving Q x , for instance in quantum correlators, will eventually depend on that hypersurface and are thus expected to break local covariance. Altogether, in the equation (7), we can approximate the local thermodynamic equilibrium operator at first order in the gradients as: where P is the total four-momentum operator and the dependence on the hypersurface of the last term was highlighted.

III. SPIN AND THERMAL SHEAR TENSOR
The mean spin polarization vector of a spin 1/2 particle can be obtained from the particle term (i.e. the future time-like part) of the Wigner function W + [23]: As pointed out in ref. [23], this formula applies to free, or quasi-free fields, therefore, in relativistic nuclear collisions, only to hadronic fields if Σ is a 3D hypersurface outside Σ F O in fig. 1. Furthermore, it is convenient to set Σ = Σ F O to calculate it most accurately. The Wigner function is the expectation value of the Wigner operator [24]: Ψ being the free Dirac field, ": :" denotes the normal ordering and θ the Heaviside step function. In the application to relativistic nuclear collisions, Ψ is to be understood as an effective hadronic field. For our purposes, the Wigner function is the expectation value of (13) with the LTE operator (6) W + ab (x, k) = Tr ρ LE W + ab (x, k) . As has been mentioned, in the hydrodynamic limit of slowly varying thermodynamic fields, one can approximate the Wigner function in x by making a Taylor expansion of the four-temperature in (6) around the point x. The contribution to the local expectation values of the gradient terms of the (11) is small compared to the contribution from the term β(x) · P , hence one can expand the exponential in (11) with the familiar techniques of linear response theory [17]: where: Thereby, the Wigner function, and the spin vector as well, in (12) will receive two linear corrections: one proportional to thermal vorticity involving correlators between the Wigner operator and the angular momentum-boost operator J x and one proportional to the thermal shear tensor, involving correlators between the Wigner operator and the operator Q x . The thermal shear tensor contribution to the spin vector has been usually neglected, for a twofold reason: first, it certainly vanishes at global equilibrium because of Killing condition and, secondly, the symmetric part of the gradients are cancelled by the Levi-Civita tensor in the formula (3). However, we will show, that a combination of thermal shear tensor and momenta eventually survives and gives rise to a contribution which can be numerically important, especially for a fluid which is not yet very close to global equilibrium.
To make calculations more compact, we will not separate the symmetric and antisymmetric part of the Taylor expansion of β like in eq. (8) and study the linear response in terms of the perturbation: where: Therefore, the particle term of the Wigner function at LTE can be approximated by: with: where with · · · β(x) we denote the thermal expectation values calculated at the homogeneous global thermodynamic equilibrium, i.e. with the density operator: The subscript c on the thermal average in (16) denotes the connected part of the correlator, that is, for the simplest case of two operators: The two terms of the right hand side of the eq. (15) can be evaluated with standard techniques (see Appendix B) and turn out to be: where T is We can now replace the ∆β in eq. (17) by using eq. (14): It is convenient to work out the integration over Σ in the eq. (19) by using the Gauss' theorem, and splitting it into an integral over a flat 3D hypersurface Σ B and a 4D integral over the region Ω B encompassed by Σ and Σ B ; for instance, for heavy ion collisions applications, it is convenient to choose Σ = Σ F O , see fig. 1. The formula (12) is thus the sum of a 4D integral and a 3D boundary term: where: and Consider the 4D integral first. Assuming that the region Ω B is large enough, we can approximate it with: where δt is the temporal extent of the region Ω B . Hence: Plugging this expression in the (20), taking into account that and keeping in mind the (18), we readily find that tr γ µ γ 5 ∆ Ω W + (x, k) = 0, so that the 4D integral does not contribute to the spin vector. We note that such a result is naturally expected for the angular momentum-boost operators, as they are a conserved charge independent of the integration hypersurface, but it is not obvious for the Q x pseudo-tensor. Indeed, the vanishing of the volume term is a specific result owing to the choice of the region where the Gauss theorem has been applied and the observable, the spin polarization vector.
We can now move on to the Σ B term. We start to evaluate the numerator: and the denominator: of the second term on the right hand side of (20). If the hypersurface Σ B is large compared to the other scales, it can be approximated with an unbounded hyperplane and so: wheret is the unit vector normal to Σ B (which corresponds to the time direction in the QGP frame of figure 1) and ∆ µν = η µν −t µtν , ∆t = (y − x) ·t and y ·t is constant in Σ B by definition. Notice that both terms in the integral of N contain the momenta p and p contracted with the Levi-Civita tensor and that the second term in (22) sets p = p in N after integrating p ; we therefore conclude that the second term in (22) does not bring any contribution. Therefore: Integrating by parts in p and taking advantage of the vanishing of the square bracket for p = p, we obtain: where we have to take into account that ∂p 0 /∂p κ is non-trivial being p on-shell. We can then integrate in p getting: where k is on-shell. Notice that the index κ is in fact a spatial index, therefore: Using the previous derivative and the decompositions we can rewrite N as: Now we can split the gradient of β as the sum of the thermal vorticity and the thermal shear: thus obtaining the linear contributions from the angular momentum-boost operators J x and the operator Q x in the (11). For the thermal vorticity we obtain: which is the known expression in eq. (3). Similarly for the thermal shear term, simple calculations yield: whence it follows, replacing the four-momentum k with p: which is, in general, non-vanishing. As has been mentioned, the integration hypersurface Σ, in relativistic nuclear collisions, is the freeze-out Σ F O (see figure 1). This term is a new, non-dissipative contribution to the spin polarization vector at local thermodynamic equilibrium to be added to the (3).

IV. DISCUSSION AND CONCLUSIONS
The striking difference between the (3) and the new term (23) is that the latter apparently breaks covariance, for the presence of thet ν , the time direction in the QGP frame (see figure 1). The reason is the inevitable dependence of Q x operator (9) on the particular 3D hypersurface of integration, as it has been discussed in section II. The appearance of a particular vectort, which is, in a sense, the best approximation of the unit vector perpendicular to the hypersurface Σ F O , is the telltale sign of this dependence.
The tensor ξ can be covariantly decomposed along the four-velocity of the fluid into: In the (24), denoting by ∇ µ = ∂ µ − u µ u · ∂ and ∆ µν = g µν − u µ u ν , A = u · ∂u is the acceleration field, σ is the properly called shear tensor: and θ = ∇ · u the expansion scalar. All of these terms can contribute to the spin vector (23), however only one of them, namely the temperature gradient, has a non-vanishing non-relativistic limit: while all remaining terms, including the spin-shear coupling, are purely relativistic. Such a term was already noticed in the non-relativistic limit of thermal vorticity [25], and provides an equal contribution. This will be the subject of further work. The additional term of spin polarization vector (23) is linear in the gradients of the thermodynamic fields and can then play a major role in driving the local spin polarization pattern in relativistic heavy ion collisions. We will show its numerical impact in a forthcoming paper [26].

ACKNOWLEDGMENTS
While we were finalizing this work, we got to know of a simultaneous similar study [27]. M.B. is supported by the Florence University fellowship Effetti quantistici nei fluidi relativistici.
Appendix A: Canonical stress-energy tensor and operator Qx The purpose of this section is to show that, for the free Dirac field, the operator Q x in eq. (9) with the Belinfante stress-energy tensor: is the same if we replaced T B with the canonical stress-energy tensor T C of the free Dirac field. Starting from the pseudo-gauge transformation relation between T B and T C : where S is the canonical spin tensor of the free Dirac field. This is known to be dual to the axial current, hence it is completely antisymmetric in the three indices; hence, the above transformation formula simplifies to: By plugging this formula into Q x , we obtain: We now focus on the integral term involving the spin tensor. Integrating by parts we get: The second integral on the right hand side vanishes because of the anti-symmetry of indices, so we are left with a 3D integral which can be turned into a surface 2D integral by means of the relativistic Stokes theorem: This integral can be made vanishing by enforcing suitable boundary conditions on the Dirac field; if Σ = Σ F O ∪ σ ± (see figure 1) and discussion in section II), usually anti-periodic boundary conditions for a compact hypersurface Σ are enforced before taking the limit to infinity. Therefore, only the first integral in eq. (A2) contributes to the operator Q x , which is the same in form as in eq. (A1).