Polarized baryon production in heavy ion collisions: an analytic hydrodynamical study

We utilize known exact analytic solutions of perfect fluid hydrodynamics to analytically calculate the polarization of baryons produced in heavy ion collisions. Assuming local thermodynamical equilibrium also for spin degrees of freedom, baryons get a net polarization at their formation (freeze-out). This polarization depends on the time evolution of the Quark-Gluon Plasma (QGP), which can be described as an almost perfect fluid. By using exact analytic solutions, we thus can analyze the necessity of rotation (and vorticity) for non-zero net polarization. In this paper we give the first analytical calculations for the polarization four-vector. We use two hydrodynamical solutions; one is the spherically symmetric Hubble flow (a somewhat oversimplified model, to demonstrate the methodology). The other solution which we use is a somewhat more involved one that corresponds to a rotating and accelerating expansion, and is thus well suited to investigate some main features of the time evolution of the QGP created in peripheral heavy-ion collisions (although there are still many numerous features of a real collision geometry that are beyond the reach of this simple model). Finally we illustrate and discuss our results on the polarization.


Introduction
Our aim is to give analytical results for the polarization four-vector of massive spin 1/2 particles produced in heavy-ion collisions, from hydrodynamical models. The motivation for this work is the recently observed non-vanishing polarization of Λ baryons at the STAR experiment [1,2] that hints at local thermal equilibrium also for spin degrees of freedom in the Quark Gluon Plasma (QGP) produced in heavy-ion collisions. The assumption of thermal equilibration for spin is at the core of the current understanding of polarization of particles produced from a thermal ensemble (such as the QGP), and almost all studies aimed at describing it in terms of collective models utilize the formula derived from this assumption by Becattini et al. [3].
Although many numerical hydrodynamical models do indeed predict non-zero polarization of produced spin 1/2 particles [4][5][6][7], a clear connection between the initial state, the final state and the observable polarization is to be expected from analytical studies, on which topic we do the first calculations here (to our best knowledge).
The observable quantities at the final state of the hydrodynamical evolution can be described by utilizing kinetic theory. At local thermodynamical equilibrium, for spin 1/2 particles such a description can be based on the the Fermi-Dirac distribution: where p µ is the four-momentum of the produced particle, and u µ (x), µ(x) and T(x) are the four-velocity, the chemical potential, and the temperature field of the fluid, respectively.

arXiv:1812.05587v2 [hep-ph] 19 Mar 2019
Assuming local thermal equlibrium for the spin degrees of freedom, for the space-time-and momentum-dependent polarization four-vector S(x, p) µ of the produced particles the following formula is given in Ref. [3]: where m is the mass of the investigated particle, and the inverse temperature field β µ =u µ /T(x) is introduced. Here ε µνρσ is the totally antisymmetric Levi-Civita-symbol; the ε 0123 = 1 convention is used. In this paper we use this formula to calculate the polarization four-vector at the freeze-out from analytical relativistic hydrodynamical solutions. The general consensus is that the appearance of polarization strongly depends on the rotation of the expanding QGP fireball. However, the Equation of State (EoS) of the QGP influences the rotation, thus by measuring the polarization, we can get information about the EoS of the QGP. Analytic hydrodynamic calculations may provide special insight by yielding analytic formulas for the connections of the aforementioned physical quantities.
We investigate two hydrodynamical solutions: the spherically symmetric Hubble flow [11,12] and a rotating and accelerating solution (first reported in Ref. [14], then in a different context in [15]). We expect to obtain zero polarization in the case of the spherical symmetric Hubble-flow as it has no rotation, so the study of this solution can be regarded as a simple cross-check of our methodology. The second one, however, being a rotating and expanding solution, could be a well usable model of peripheral heavy-ion collisions, and it is expected that one gets non-zero polarization out of it. Thus this rotating expanding solution constitutes the core point of the reported work.

Basic equations and assumptions
We use the c=1 notation. Let us denote the space-time coordinate by x µ ≡(t, r), and the Minkowskian metric tensor by g µν =diag(1, −1, −1, −1). The convention for the Levi-Civita symbol is ε 0123 =1. Greek letters denote Lorentz indices, Latin letters denote three-vector indices. For repeated Greek indices we use the Einstein summation convention. We denote the space dimension by d; this implies g µ µ = d+1. In reality, d = 3, but it is useful to retain the d notation wherever possible, in order to see if the reason for a specific numeric constant in the formulas is the dimensionality of space. The four-velocity of the fluid is u µ =γ(1, v), where γ= √ 1−v 2 is the Lorentz factor. The velocity three-vector is then v=u k /u 0 . With p µ we denote the four-momentum of a produced particle; we also use the three-momentum p, whose magnitude we simply denote by p (whenever there is no risk of confusion). The energy of the particle is denoted by E; the mass shell condition then reads as E = p 2 +m 2 , with m being the particle mass.
The usability of hydrodynamics in heavy ion physics phenomenology relies on the assumption of local thermodynamical equilibrium of the matter. For describing particles with spin 1/2, we use the source function as written up in Eq. (1). Hadronic final state observables can be then calculated by integrating over the freeze-out hypersurface; e.g. in the case of the invariant momentum distribution, the driving formula is Here d 3 Σ ν is the 3-dimensional vectorial integration measure of the freeze-out hypersurface; the appearance of which is the so-called Cooper-Frye prescription [8] for calculating the invariant momentum distribution. Of the two solutions (mentioned above) which we investigate in this work, in the case of the rotating and expanding accelerating solution, we also calculate the invariant momentum distribution, as this has not been done before. The formula given in Ref. [3] for the polarization of spin 1/2 particles, as written up in Eq. (2), may be utilized for any given β µ =u µ /T field that one gets from a given solution of the hydrodynamical equations. We are interested in calculating the polarization at the final state of the hydrodynamical evolution, so we must integrate the S(x, p) µ field over the freeze-out hypersurface. The formula to be analyzed further, that is, that for the observed polarization S(p) µ of particles with momentum p, thus becomes as written up e.g. in [7]. For being able to perform analytical calculations, we have to make some assumptions. We use saddle point integration, in which one assumes that the integrand is of the form f (r)g(r), where f (r) is a slowly changing function, while g(r) has a unique and sharp maximum; then the integral can be calculated with a Gaussian approximation as that is, R 0 is the location of the unique maximum of g(r) and M is the second derivative matrix.
Another assumption concerns the expression of S(x, p) µ , Eq. (2): if the exponent in the Fermi-Dirac distribution is large (i.e. phase space occupancy is small), we can use the Maxwell-Boltzmann distribution instead: Here g is the spin-degeneracy factor; for spin 1/2 baryons, g = 2.
In high energy heavy ion phenomenology (when the collision energy is high enough, say for collisions at RHIC or LHC), the µ/T factor can (and usually is) neglected; we use this approximation here 1 . With this we have If the Maxwell-Boltzmann approximation is justified, it means that f (x, p) 1 indeed, and then also Eqs. (2) and (4) become simpler: and in the saddle-point approximation, the polarization of particles with momentum p becomes simply since in the saddle-point approximation, in the numerator of Eq. (4), S(x, p) µ can be considered the "smooth" function, and the determinant factors cancel. 1 The vanishing of µ can also be interpreted as an absence of a conserved particle number density n. All our conclusions would change only by a proportionality factor if we said µ/T = const instead of µ/T = 0; if µ = 0, we would have had to introduce n. Depending on the equation of state of the matter (one that also contains the conserved particle density n), one could write the f (x, p) function in another form, where the normalization dp f (x, p) = n(x) is evident. For example, if one chooses an ultra-relativistic ideal gas, with p = nT, ε = κ p, with κ = d as EoS, one has g (2πh) d e µ/T = n 4πT 3 . Indeed, in the solutions discussed below, µ/T =const is satisfied, which means n ∝ T d , which is the well-known condition for an adiabatic expansion.

Some exact hydrodynamical solutions and polarization
In this section we first specify and recapitulate the investigated hydrodynamical solutions, then give the analytical formulas for the polarization four-vector calculated from them. The equations of perfect fluid relativistic hydrodynamics utilized here are n∂ µ u µ = −u µ ∂ µ n (particle number/charge conservation), and we specify the simple ε = κ p equation of state here. (The notations: ε, p and n are the energy density, pressure and particle number density, respectively.) Concerning the n density: if it is assumed to be non-vanishing, we set the EoS as p = nT. However, the solutions presented below are valid also if n = 0 (ie. if µ = 0). So the expressions for n that we recapitulate for the solutions can be regarded as supplemental to the solutions that work for µ = 0. We also note that there is recent development on taking the effect that polarization of the constituents of the fluid has on the fluid dynamics itself [9], along with some numerical calculations of how this modified hydrodynamical picture affects final state polarization [10]. We do not investigate this possibility here; we restrict ourselves to the simple and well-known basic equations witten up above.

Hubble flow
We do not go into the details about the method to find or verify that the solutions presented below are indeed solutions of the perfect fluid hydrodynamical equations; we refer back to the original publications of the solutions.
We investigate the Hubble-like relativistic hydrodynamical solution first fully described in Ref. [11]. This solution has the following velocity, particle density and temperature fields: where τ = √ t 2 − r 2 , and κ is the inverse square speed of sound (constant in the case of this exact solution). The κ=3 case corresponds to ultrarelativistic ideal gas, κ=3/2 corresponds to a non-relativistic gas; however, this solution is valid for any arbitrary constant κ value 2 .
To calculate the polarization four-vector, as of now we investigate the simplest case, the spherical symmetric expansion. For the freeze-out hypersurface the τ=τ 0 =const. hypersurface is chosen (which in the case of the investigated solution equals the constant temperature freeze-out hypersurface), and a given point of this hypersurface can be parametrized simply by the r coordinate three-vector, and 2 We note that a more general class of solutions is possible [11][12][13] in which the temperature and density fields are supplemented with an arbitrary V function of a "scaling variable" S: and the S variable is any function of S x , S y , and S z : HereẊ 0 ,Ẏ 0 andŻ 0 are arbitrary constants. In the given example, the S = const surfaces are ellipsoids, andẊ 0 ,Ẏ 0 ,Ż 0 are time derivatives of the principal axes of them the time coordinate on the hypersurface is t(r)≡ τ 2 0 +r 2 . The integration measure and the resulting expression for the Cooper-Frye formula can then be written as As we are discussing massive particles, this integral always exists. The T 0 constant (an arbitrary parameter of the solution) can be taken simply as the temperature at freeze-out; we did so. The position of the saddle-point (R 0 ) as well as the second derivative matrix M kl is calculated as With this we can get an approximation for the invariant single-particle momentum distribution: The formula is independent of momentum. This was expected because this hydrodynamical solution (in the V (S)=1 case) is boost invariant. To use (9) to determine the polarization four-vector in the hydrodynamical solution of the Hubble-flow, first we give the expression for the ∂ ν β ρ derivative: Then for the time component we get: as ε 0ikl is antisymmetric whereas g ik and r i r k are symmetric to the change in the i ↔ k indices. Similarly for the spatial coordinates: In conclusion, the polarization four-vector in the spherical symmetric Hubble-flow is which is consistent with our expectations.

Rotating and accelerating expanding solution
Another hydrodynamical solution of particular interest to us is a rotating and accelerating expanding solution, first written up in Ref. [14]. This solution has the following velocity, temperature and particle density profiles: where ρ 0 and τ 0 are arbitrary parameters and Ω is an arbitrary angular velocity three-vector that indicates the axis and magnitude of rotation. The ρ 0 parameter tells about the initial spatial extent of the expanding matter, however, the τ 0 parameter is just there for the sake of consistency of physical units; in this way, the unit of Ω is c/fm, as it should be for an angular velocity-like quantity 3 , and T 0 is a temperature constant. In the case of Ω=0, we recover an acceleratingly expanding but non-rotating spherically symmetric solution.
It is convenient to write up this solution with the following notation: To calculate final state observables, we choose the constant proper time (τ 0 =const) hypersurface here as well. The solution itself allows for a re-scaling of the arbitrary constants in the formulas; just as in the previous case, here too we can treat the T 0 quantity as the temperature at freeze-out (at the r = 0 center of the expanding matter). We use the notation introduced in Eq. (11) for the Maxwell-Boltzmann distribution. To derive the saddle point for the calculation of the polarization four-vector, we shall use the expression of the invariant momentum spectrum: This integral always exists (in the case of massive particles). In order to utilize the saddle-point integration method, we determine the position of the saddle point (R 0 ) and the second derivative matrix at the saddle point: We leave the detailed calculations to Appendix .1; the results are the following. The R 0 saddle point (for a given p momentum) is in the plane spanned by the p and p×Ω vectors. In the following we use thep ≡ p/p notation for the unit vector pointing in the direction of p. For the saddle point we get 3 Here we changed the notation of Ref. [14]. The rather unfortunate B notation used there is now written as τ 2 0 Ω.
Concerning the second derivative matrix, we need it only for the calculation of the invariant momentum distribution, where its determinant is invoked. It turns out that this quantity is Using this result, we get the invariant single-particle momentum distribution 4 as Equivalently, by defining a "local slope" T eff , the result can be expressed as Proceeding to the polarization of the produced baryons, we calculate the derivative of the inverse temperature field for this solution from the form given in Eq. (19), then substitute it into the expression of the polarization, Eq. (9). The result is (The second term was cancelled owing to the symmetry of g νρ and the antisymmetry of ε µνρσ , and x µ is understood as the four-coordinate of the freeze-out hypersurface whose three-coordinate is the r = R 0 three-vector). Remembering the expression of the introduced F µν tensor and b µ vector from Eq. (19), in particular that F 0k = 0, and b k = 0, we get the following expressions for the the time-like and space-like components: Summarizing this result, the polarization four-vector for the investigated rotating and accelerating expanding solution is the following: In the case of Ω = 0, there is no rotation, and we get S(p) µ =0. In this model thus polarization is very transparently connected to the presence of rotation. 4 This has not yet been calculated for this hydrodynamical solution.
It is useful to transform the polarization four-vector into the rest frame of the particle. The result is 5 , with (r.f. standing for "rest frame"): We can also compute the helicity of the produced spin 1/2 particles in this solution from this formula (the S polarization vector is taken in the laboratory frame):

Illustration and discussion
In this section, we would like to illustrate our simple analytic results for the polarization vector. We use the same type of plots that was used to visualize some existing numerical simulations (e.g. those presented in Ref. [7]). We plot the components of the polarization vector with respect to the momentum components in the transverse plane (that is, w.r.t. p x and p y ). On Fig. 1 we plot the polarization vector in the laboratory frame. For the sake of plotting, the mass of the Λ baryon (m Λ c 2 = 1115 MeV) was chosen. For the sake of this illustration, we chose a moderate value for the magnitude of the Ω vector as |Ω| = 0.1 c/fm. In our case, as a special coincidence owing purely to the specific algebraic form of the presented analytic solution, it turned out that the polarization in the rest frame of the produced baryons is independent of momentum p; see Eq. (32). This coincidence is expected to be relieved in the case of more involved (complicated) solutions (that are left for future investigations). Fig. 2 nevertheless shows the value of the S y component in the baryon rest frame. The helicity of the produced baryons (being proportional to the pS scalar product), however, does depend on the momentum, even in the case of our very simple solution. We plot it on Fig. 3; with the same parameter values as in the foregoing two plots.

Summary and outlook
In this paper we gave the first analytical formulas for the polarization of baryons produced from a thermal ensemble corresponding to rotating and expanding exact hydrodynamical solutions. These arise as descriptions of the final state of non-central high energy heavy-ion collisions. We investigated two exact relativistic hydrodynamical solutions. One was the spherically symmetric Hubble flow (an overly simplistic one, the study of which can be regarded as a check of the methodology), in which the polarization turns out to be exactly zero (as it is naturally expected from symmetry considerations). The other solution we investigated is a one describing rotating and accelerating expansion. In this case we got the first ever analytical formulas that connect dynamical quantities of the expansion (i.e. magnitude of rotation, acceleration, etc) with the observable final state polarization of spin 1/2 particles (baryons), which turns out to be non-zero in this case.
Our results are simple and straightforward. Nevertheless, many more solutions (more involved ones) as well as more complicated final state parametrizations can be investigated in the future. The calculations presented here yield the first results in terms of exact formulas for the polarization; more refined future studies are needed to disentangle the effects that rotation, acceleration and temperature 5 The Lorentz matrix performing this boost transformation is the following (in usual 1+3 dimensional block matrix notation): where E and p could be parametrized with the velocity parameter χ as E = m cosh χ and p = m sinh χ, respectively. It indeed can be checked that this matrix takes the (E, p) four-momentum vector into (m, 0), as it should.   where we temporarily introduced the A≡ τ 2 0 +α 2 p 2 + β 2 τ 4 0 (p 2 Ω 2 −(pΩ) 2 ) notation. Because of the orthogonality of p and p×Ω, both sides here have to vanish identically, from which we get One divides these equations to obtain a simple relation, the substituting back one gets a quadratic equation for β, the solution of which is where we used the E 2 =p 2 +m 2 relation. To find α we substitute this back into the expression of A: Using the above expression of β (with the yet undetermined sign) we get 1−4p 2 β 2 = − 2m p 2 (m±E), and see that the expression for α will be valid only in the case when 1−4β 2 p 2 > 0, thus conclude that the bottom sign is the proper choice. We thus arrive at the following expressions: where we use the notation A as above. We should use the expression of R 0 as calculated above. The determinant of this M matrix is the product of its eigenvalues. In our case the particular spatial directions are: p, p×Ω, and the vector orthogonal to both these, that is, p×(p×Ω). One recognizes that the vector p×(p×Ω) is an eigenvector of the M second derivative matrix: M p×(p×Ω) = · · · = 1 β p×(p×Ω).
The corresponding eigenvalue is thus 1/β. Owing to the symmetric nature of M, the other two eigenvectors must be in the orthogonal complementer subspace of this vector, so they are linear combinations of p and p×Ω. Let us thus look for these eigenvectors in the form a=µp+νr, with yet to be determined µ and ν coefficients. Substituting this expression, we get where λ is the eigenvalue (the values of which we are looking for). By substituting the assumed form of a and inferring the components of this equation in the p and p×Ω directions, we get the following equation for the µ and ν coefficients: We immediately infer the product of the two λ 1,2 eigenvalues as the determinant of this 2×2 matrix. Taking the third eigenvalue (calculated above) into account, after some simplifications, we indeed get the following expression for the determinant of the M matrix (the expression we used in Eq. (25)): det M = 1 T 0 τ 2 0 3 32pm 2 (m+ p 2 +m 2 ). (A9)