Study of (cid:2) polarization in relativistic nuclear collisions at √ s NN = 7 . 7–200 GeV

We present a calculation of the global polarization of (cid:2) hyperons in relativistic Au–Au collisions at RHIC Beam Energy Scan range √ s NN = 7 . 7 − 200 GeV with a 3 + 1-dimensional cascade + viscous hydro model, UrQMD + vHLLE. Within this model, the mean polarization of (cid:2) in the out-of-plane direction is predicted to decrease rapidly with collision energy from a top value of about 2% at the lowest energy examined. We explore the connection between the polarization signal and thermal vorticity and estimate the feed-down contribution to (cid:2) polarization due to the decay of higher mass hyperons.


Introduction
Particles produced in relativistic heavy ion collisions are expected to be polarized in peripheral collisions because of angular momentum conservation. At finite impact parameter, the Quark-Gluon Plasma (QGP) has a finite angular momentum perpendicular to the reaction plane and some fraction thereof may be converted into spin of final state hadrons. Therefore, measured particles may show a finite mean global polarization along the angular momentum direction.
Early estimates of this effect [1] were based on the general idea that polarized quarks in the QGP stage of the production process would eventually give rise to polarized hadrons, making it possible to predict qualitative features of the final hadrons polarization. It was then proposed [2,3] that polarization can be calculated assuming that the spin degrees of freedom are at local thermodynamical equilibrium at the hadronization stage in much the same way as the momentum degrees of freedom. In other words, polarization can be predicted finding the appropriate extension of the familiar Cooper-Frye formula to particles with spin. A specific a e-mail: yu.karpenko@gmail.com derivation was presented in Refs. [2,4], where it was pointed out that the hydrodynamical quantity steering the polarization is the thermal vorticity, that is, (minus) the antisymmetric part of the gradient of the four-temperature field β = (1/T )u where T is the proper temperature and u the hydrodynamic velocity: Particularly, the first-order expansion of the polarization in terms of thermal vorticity was obtained in Ref. [4] for hadrons with spin 1/2 (lately recovered with a different method in Ref. [5]), yet its extension to higher spins could be derived from the corresponding global equilibrium expression [2]. This theoretical work made it possible to make definite quantitative predictions of global polarization in nuclear collisions from hydrodynamic calculations, with a resulting mean value ranging from some permille to some percent [6][7][8], with an apparently strong dependence on the initial conditions, particularly on the initial longitudinal velocity field. Calculations of vorticity in relativistic heavy ion collisionswhich could then be turned into a polarization map -were also recently presented in Ref. [9]. To complete the theoretical overview on the subject, it should be pointed out that different approaches, as well as additional mechanisms, to the polarization in relativistic nuclear collisions were proposed in Refs. [10][11][12][13][14].
From the experimental viewpoint, while the early measurement by the STAR experiment in Au-Au collisions at √ s NN = 200 GeV [15] only set an upper limit on the polarization of 0.02, a first strong evidence of non-vanishing global polarization showing the predicted features within the hydrodynamical model (orthogonal to the reaction plane, equal sign for particle and antiparticle) has been recently reported by the STAR experiment [16,17] at the lower ener-gies of the RHIC Beam Energy Scan (BES), with a top value at the lowest energy point √ s NN = 7.7 GeV of some percent. This finding certainly opens a new direction in the physics of the Quark-Gluon Plasma and relativistic heavy ion collisions and, in the short run, calls for new numerical calculations of thermal vorticity and polarization in the BES energy range, what is the primary purpose of this work. As has been mentioned, previous quantitative estimates of polarization in hydrodynamic models vary from some permille to some percent, with strong dependence on the collision model and, specifically, on the hydrodynamic initial conditions. In [6], simulation of noncentral √ s NN = 200 GeV Au-Au collisions at RHIC with initial state from Yang-Mills dynamics followed by a 3-dimensional ideal hydro expansion implying an initial non-vanishing space vorticity vector [18] resulted in polarization of lowp T of few percent magnitude, parallel to the direction of total angular momentum of the fireball, whereas for few GeV p T it reaches about 8%. In Ref. [7], at the same collision energy, a 3-dimensional hydrodynamic expansion with initial state from the optical Glauber model with parametrized space-time rapidity dependence chosen to reproduce the momentum rapidity distributions of hadrons and directed flow results [19] resulted in much lower values of polarization: about 0.2% for lowp T and up to 1.5% for the longitudinal (along the beam) component of polarization of highp T . In a recent work [8] an event-by-event 3-dimensional viscous hydrodynamics with initial state from AMPT model results in a similar few per-mille average global polarization for A + A collisions at √ s NN = 62.4, 200 and 2760 GeV and the correlation of polarizations of pairs is studied.
Here, we present the simulations of Au-Au collisions at the RHIC BES energies, √ s NN = 7.7 . . . 200 GeV using a 3-dimensional viscous hydrodynamic model vHLLE with initial state from the UrQMD model. Such a hybrid approach coupled to final state UrQMD has been tuned to reproduce the basic hadron observables in relativistic heavy ion collisions, that is, (pseudo)rapidity and transverse momentum distributions and elliptic flow coefficients. We calculate the polarization of hyperons produced out of the fluid (socalled thermal ). Besides, we estimate, for the first time, the contribution to the polarization stemming from the decays of likewise polarized 0 , (1385), as well as other and resonances up to (1670). The paper is organized as follows: in Sect. 2 we briefly describe the model used to simulate nuclear collisions in the examined energy range; in Sect. 3 we summarize the main definitions concerning spin and polarization in the relativistic regime, as well as the used formulas to calculate the polarization; in Sect. 4 we present the results for the polarization of baryons in Au-Au collisions at BES energies and discuss their interpretation in terms of vorticity patterns in the hydrodynamic expansion; in Sect. 5 we estimate the correc-tions due to the resonance decays. Conclusions are drawn in Sect. 6.

Nuclear collision model description
The full cascade + viscous hydro + cascade model used in our studies is described in detail in Ref. [20], herein we only summarize its main features. At lower collision energies the colliding nuclei do not look like thin "pancakes" because of weaker Lorentz contraction; also, partonic models of the initial state (CGC, IP-Glasma) gradually lose their applicability in this regime. The longitudinal boost invariance is not a good approximation anymore, therefore one needs to simulate a 3-dimensional hydrodynamic expansion. Such considerations motivate us to choose UrQMD model to describe the dynamics of the initial state from the first nucleon-nucleon collisions until a hypersurface of constant Bjorken proper time τ = √ t 2 − z 2 = τ 0 to provide a 3 dimensional initial state for subsequent hydrodynamic stage. 1 For most of the collision energies, the value of τ 0 is chosen to correspond to an average time when the two incoming nuclei have passed through each other, τ 0 = 2R/(γ v z ).
At the τ = τ 0 hypersurface the transition to a fluid description occurs: energies and momenta of particles crossing the hypersurface are distributed to the hydrodynamic cells around the positions of particles according to a Gaussian profile. The contribution of each particle to a fluid cell {i, j, k} is given by: where x i , y j , η k are the differences between particle's position and coordinates of the center of the hydrodynamic cell {i, j, k} in the transverse plane and space-time rapidity, R ⊥ and R η are the coarse-graining (smearing) parameters and γ η = cosh(y p − η) is the longitudinal Lorentz factor of the particle as seen in a frame moving with the rapidity η; the normalization constant C is calculated in order to conserve energy/momentum in this transition. With such initial conditions the 3-dimensional viscous hydrodynamic evolution in the Israel-Stewart formulation starts, which is numerically solved with vHLLE code [21]. In particular, we work in a Landau frame, and we neglect baryon and electric charge diffusion currents. Whereas the bulk viscosity is set to zero, ζ /s = 0, for the shear-stress tensor π μν the following evolution equation is solved: where the relaxation time of the shear-stress tensor is set to τ π = 5η/(T s), the brackets denote the traceless part and the part orthogonal to u μ of the tensor and π μν NS is the Navier-Stokes value of the shear-stress tensor.
Finally, we set the transition from fluid to particle description (a.k.a. particlization) to happen at a certain energy density = sw . The elements of particlization hypersurface are computed in the course of hydrodynamic evolution by means of the CORNELIUS subroutine [22]. To calculate basic hadronic observables, the Monte Carlo hadron sampling at the particlization surface is performed using basic Cooper-Frye formula with standard quadratic ansatz for the shear viscous corrections. The sampled hadrons are fed into UrQMD cascade to calculate the final state hadronic interactions.
It has been shown in Ref. [20] that in such a model it is not possible to reproduce the basic bulk hadronic observables -(pseudo)rapidity distributions, transverse momentum spectra and elliptic flow coefficients -simultaneously in the collision energy range √ s NN = 7.7 . . . 200 GeV if the parameters of the model (except the hydro starting time τ 0 ) are fixed at constant values. However, a reasonable reproduction of the experimental data has been achieved when the parameters of the model were chosen to depend monotonically on the collision energy as it is shown in Table 1. This was obtained when the particlization energy density was fixed at sw = 0.5 GeV/fm 3 for the whole collision energy range.
The UrQMD cascade does not treat polarization of hadrons, therefore we calculate the polarization of particles produced at the particlization surface only. To to so, have replaced the Monte Carlo hadron sampling with a direct calculation based on the Eq. (10) applied on particlization surfaces from event-by-event hydrodynamics. The final state hadronic cascade (UrQMD) is not used in the present study.

Spin and polarization
We summarize in this section the basic definitions concerning massive relativistic particles with spin.
The starting point is the so-called Pauli-Lubanski vector operator W μ , which is defined as follows: where J and P are the angular momentum-boost operators and four-momentum operators for a single particle. The proper spin four-vector operator is simply the adimensional vector obtained by dividing W by the mass, that is, It can easily be shown [23] that According to the latter relation, the spin vector is spacelike and has only three independent components. Therefore, for single particle states with definite four-momentum p it can be decomposed [24] along three spacelike vectors n i ( p) with i = 1, 2, 3 orthogonal to p: It can be shown that the operators S i ( p) with i = 1, 2, 3 obey the well known SU(2) commutation relations and they are indeed the generators of the little group, the group of transformations leaving p invariant for a massive particle. In other, maybe simpler, words they correspond to the familiar spin operators of the non-relativistic quantum mechanics in the particle's rest frame. Also, the S μ S μ operator is a Casimir of the full Poincaré group whose eigenvalue is S(S+1) where S is the spin of the particle. We denote the mean spin vector the four-vector obtained by taking both the quantum and the thermal average of S, that is, and each component takes on values in the interval (−S, S). Indeed, the properly called polarization vector is obtained by normalizing μ to the spin of the particle, that is, hence each component takes on values in the interval (−1, 1). For a multi-particle system, the calculation of the polarization of particles with four-momentum p requires the decomposition of the mean total angular momentum into momentum modes. In general, this requires the knowledge of the Wigner function, which allows one to express the mean values of operators as integrals over space-time and fourmomentum space. By using such theoretical tools, the mean spin vector of spin 1/2 particles with four-momentum p, produced around point x on particlization hypersurface, at the leading order in the thermal vorticity reads [4] a formula recovered in Ref. [5]. In Eq. (9) f (x, p) is the Fermi-Dirac distribution and is given by Eq.
(1) at the point x.
In the hydrodynamic picture of heavy ion collisions, particles with a given momentum are produced across the entire particlization hypersurface. Therefore to calculate the relativistic mean spin vector of a given particle species with given momentum, one has to integrate the above expression over the particlization hypersurface [4]: The mean (i.e. momentum average) spin vector of all particles of given species can be expressed as where N = d 3 p , p) is the average number of particles produced at the particlization surface.
In the experiment, the polarization is measured in its rest frame, therefore one can derive the expression for the mean polarization vector in the rest frame from (11) taking into account Lorentz invariance of most of the terms in it: where an asterisk denotes a quantity in the rest frame of particle.
As has been mentioned, Eq. (9) applies to spin 1/2 particles. However, a very plausible extension to higher spins can be obtained by taking into account that the expression of the spin vector for massive particle of any spin in the Boltzmann statistics is known at global equilibrium with rotation [2]; this is discussed in detail in Ref. [25]. The extension of Eq. (9), in the limit of Boltzmann statistics, reads

Coordinate system
We use the following notation for the reference frame in the simulations: the z axis points in the direction of the beam, and the x axis is parallel to the impact parameter. Thus xz is the reaction plane, and y axis points out of the reaction plane, oppositely oriented to the total angular momentum J. For the components of the polarization vector (8) we use physically motivated aliases: P || is component in the direction of beam (z), P b is component in direction of impact parameter (x) and P J is component in direction of the total angular momentum of the fireball J (i.e. −y) (Fig. 1).

Polarization of thermal as a function of transverse momentum
We started simulating semi-central Au-Au collisions at one particular collision energy in the middle of the Beam Energy Scan range, √ s NN = 19.6 GeV. We ran 1000 event-byevent hydro calculations with initial state corresponding to 40-50% centrality class (impact parameter range b = 9.3 − −10.4 fm) and calculated event-averaged denominator and numerator of Eq. (10) in the p x p y plane at zero momentum space rapidity. The resulting components of momentum differential polarization vector [see Eq. (8) for a definition] are shown on Fig. 2. One can see that the p T dependence of the polarization has a pattern which is similar to the one obtained in a 3 + 1D hydrodynamic simulation of √ s NN = 200 GeV Au-Au system with ECHO-QGP code [7,26], including the signs of polarization in different quadrants of the p x p y plane. On this plot, we extend the calculation up to p ⊥ ≈5.5 GeV, which is well beyond the applicability of hydrodynamics, to show the trends in the momentum dependence. production decreases exponentially with p ⊥ , therefore larger polarization of highp ⊥ does not influence the momentum-averaged polarization.  The polarization patterns in the p x p y plane reflect the corresponding patterns of the components of thermal vorticity over the particlization hypersurface. In particular, we found that the leading contribution to P b stems from the term tz p y in Eq. (9). In turn, tz , shown in left panel of Fig. 3 is a result of the interplay of ∂ t β z (acceleration of longitudinal flow and temporal gradients of temperature -conduction) and ∂ z β t (convection and conduction), according to Eq. (1). The P J component has a leading contribution from the term xz p 0 (which is also the only non-vanishing contribution at p T = 0), and xz has a rather uniform profile over the zero space-time rapidity slice of the particlization hypersurface, and the leading contribution to it comes from ∂ x u z (shear flow in z direction).
At large transverse momenta, the P || component of polarization has the largest value; however, since P b and P || flip signs in different quadrants in p x p y plane, their p T integrated values vanish, and the only nonzero component remaining is P J , which is parallel to the direction of total angular momentum J of the fireball. This is in line with the physical expectation that the global polarization has to be collinear with the vector of the total angular momentum of the system.

Collision energy dependence of the polarization
Next we ran the simulations for the full BES collision energy range. To follow recent STAR measurements, we choose 20-50% centrality bin by correspondingly chosen range of impact parameters for the initial state UrQMD calculation. The resulting p T integrated polarization is shown in Fig. 4. We observe that the mean polarization component along J , that is, P J decreases by about one order of magnitude as collision energy increases from √ s NN = 7.7 GeV to full RHIC energy, where it turns out to be consistent with the results of [7].
The fall of the out-of-plane component P J is not directly related to a change in the out-of-plane component of total angular momentum of the fireball. In fact, the total angular momentum increases as the collision energy increases, which can be seen on top panel of Fig. 5. However, the total angular momentum is not an intensive quantity like polarization; therefore, to have a better benchmark we took the ratio between the total angular momentum and the total energy, J/E, which is shown in the bottom panel of the same figure. Yet, one can see that the J/E shows only a mild decrease as the collision energy increases.
In Fig. 6 we show the distribution of the average polarization of as a function of centrality (i.e. N part ), where each point corresponds to a hydrodynamic evolution with a given fluctuating initial condition characterized by N part ; in the bottom panel one can see the corresponding distribution of total angular momentum J . We observe that the total angular momentum distribution has a maximum at certain range of N part , and drops to zero for the most central events (where the impact parameter is zero) and most peripheral ones (where the system becomes small). In contrast to that, the polarization shows a steadily increasing trend toward peripheral collisions, where it starts to fluctuate largely from event to event because of smallness of the fireball, a situation where the initial state fluctuations start to dominate in the hydrodynamic stage. As has been mentioned above, the parameters of the model are taken to monotonically depend on collision energy in order to approach the experimental data for basic hadronic observables. The question may arise whether the collision energy dependence of P J is the result of an interplay of collision energy dependencies of the parameters. We argue that it is not the case: in Fig. 7 one can see how the p T integrated polarization component P J varies at two selected collision energies, √ s NN = 7.7 and 62.4 GeV, when the granularity of the initial state controlled by R ⊥ , R η parameters, shear viscosity to entropy ratio of the fluid medium η/s and particlization energy density sw change. It turns out that a variation of R ⊥ within ±40% changes P J by ±20%, and a variation of R η by ±40% changes P J by ±25% at √ s NN = 62.4 GeV only. The variations of the remaining parameters affect P J much less. We thus conclude that the observed trend in p T integrated polarization is robust with respect to variations of parameters of the model. Now we have to understand the excitation function of the p T integrated P J , which is calculated in the model. As has been mentioned, P J at low momentum (which contributes most to the p T integrated polarization) has a dominant contribution proportional to xz p 0 . It turns out that the pattern and magnitude of xz over the particlization hypersurface change with collision energy.
We demonstrate this in Fig. 8 for two selected collision energies. For this purpose we ran two single hydrodynamic calculations with averaged initial conditions from 100 initial UrQMD simulations each. At √ s NN = 62.4 GeV, because of baryon transparency effect, the x, z components of fourtemperature vector around zero space-time rapidity are small and do not have a regular pattern, therefore the distribution of xz in the hydrodynamic cells close to particlization energy density includes both positive and negative parts, as is seen on the corresponding plot in the right column. At √ s NN = 7.7 GeV, baryon stopping results in a shear flow structure, which leads to the same (positive) sign of the xz .
In the right column of Fig. 8, we plot the corresponding xz distributions over the particlization hypersurfaces projected on the proper time axis. Generally speaking, hydrodynamic evolution tends to dilute the initial vorticities. One can see that longer hydrodynamic evolution at √ s NN = 62.4 GeV in combination with a smaller absolute value of the average initial vorticity results in an average absolute vorticity a factor 4-5 smaller at late times for √ s NN = 62.4 GeV than for √ s NN = 7.7 GeV. This results in a corresponding difference in the momentum-integrated polarization at these two energies, which is mostly determined by lowp T , which are preferentially produced from the Cooper-Frye hypersurface at late times.

Feed-down and post-hadronization interactions
Thus far, we have calculated the polarization of the produced from the plasma at the particlization stage -henceforth denoted as primary -when the fluid decouples at or right after hadronization. However, a sizable amount of the final are products of resonance decays. Indeed, as long as one is interested in the mean, momentum-integrated, spin vector in the rest frame, it can be shown that a simple linear rule applies [25], that is, where D is the daughter particle, X the parent and C a coefficient whose expression may or may not depend on the dynamical decay amplitudes. If the coefficient C does not depend on the dynamical decay amplitudes, it takes on rational values  Table 2 Polarization transfer coefficients C [see Eq. (14)] to the or hyperon (the 1/2 + state) for various strong/electromagnetic decays depending on Clebsch-Gordan coefficients, the initial values of spin and parity [25]. The values which are relevant for our calculation in various strong/electromagnetic decays with a or a hyperon in the final state are reported in Table 2; for the full derivation of the C coefficients, see Ref. [25].
A large fraction of secondary come from the strong (1385) → π and the electromagnetic 0 → γ decays 2 . We found that -in our code -the fractions of pri- 2 We denote (1385) below as * for brevity. mary , from * decays and from the decays of primary 0 are, respectively, 28, 32 and 17%, with a negligible dependence on the collision energy. This is very close to the fractions extracted from a recent analysis [27] within the statistical hadronization model: 25, 36 and 17%. The remaining 23% of the consist of multiple smaller contributions from decays of heavier resonances, the largest of which are (1405), (1520), (1600), (1660) and (1670). Some of these resonances produce in cascade decays, for example (1405) → 0 π, 0 → γ .
We start with the contribution from * , which is a J π = 3/2 + state. In this case the factor C in Eq. (14) is 1/3 (see Table 2) and, by using Eq. (13) with S = 3/2, we see that the mean spin vector of primary * is 5 times the one of the primary . Thus, the mean spin vector of from * decay is Similarly, for the 0 , which is a 1/2 + state, the coefficient C is −1/3 (see Table 2) and as the primary 0 are expected to have the same polarization as . More generally, the law (14) makes it possible to calculate the mean spin vector inherited by the hyperons also in multi-step two-body decays with a simple linear propagation rule. We have thus computed the mean spin vector of primary + feed-down as follows: where the sum goes over resonances X decaying into or 0 and a pion, N X are the primary multiplicities of resonances, C X → and C X → 0 are polarization transfer coefficients, b X → and b X → 0 are the branching ratios for decay channels yielding in and 0 , respectively. In Eq. (15) we have used the fact that nearly all 0 decay to γ with polarization transfer −1/3, which allows one to treat cascade decay contributions X → 0 → . The results are shown in Fig. 9, where the thin dotted line corresponds to feeddown contributions X = 0 , (1385) only. Surprisingly, in this case the interplay of hadron chemistry and polarization transfer in the decays results in a correction factor which varies between 0.94-0.98 in the whole collision energy range. When we take all aforementioned resonances into account (X = 0 , (1385), (1405), (1520), (1600), (1660),  Table 2, we obtain the dashed line in Fig. 9, corresponding to a 15% suppression of the mean polarization of with respect to the primary polarization. This decrease is mostly due to the increase of the denominator in Eq. (15) from the heavier resonance contributions, whereas their contributions to the numerator have opposite signs because of alternating signs of the polarization transfer coefficients.
There is, however, a further correction which is much harder to assess, i.e. post-hadronization interactions. In fact, hadronic elastic interaction may involve a spin flip which, presumably, will randomize the spin direction of primary as well as secondary particles, thus decreasing the estimated mean global polarization in Fig. 9. Indeed, in UrQMD cascade which is used to treat interactions after particlization, the cross sections of and 0 with most abundant mesons and baryons -calculated with the Additive Quark Model -are comparable to those of nucleons [28]. This implies that the do rescatter in the hadronic phase, and indeed we observe from the full cascade + hydro + cascade calculation that in the RHIC BES range only 10-15% of the primary escape the system with no further interactions 3 , until they decay into pion and proton far away from the fireball. For the present, we are not able to provide a quantitative evaluation of the rescattering effect on polarization, which assessment is left to future studies. The only safe statement for the time being is that the dashed line in Fig. 9 is an upper bound for the predicted mean global polarization within the hydrodynamical model with the specific initial conditions used in our calculation.

Conclusions
In summary, we have calculated the global polarization of hyperons produced at midrapidity in Au-Au collisions at RHIC Beam Energy Scan collision energies, √ s NN = 7.7-200 GeV, in the framework of 3-dimensional event-byevent viscous hydrodynamic model (UrQMD+vHLLE). The in-plane components of the polarization vector as a function of transverse momentum are found to have a quadrupole structure (similar to the one obtained in [7]) and can be as large as several percents for large transverse momentum. The mean, momentum-integrated polarization vector is directed parallel to the angular momentum of the fireball and its magnitude substantially increases from 0.2 to 1.8% as collision energy decreases from full RHIC energy down to √ s NN = 7.7 GeV. Such an increase is related to (1) the emerging shear flow pattern in beam direction at lower colli-sion energies related to baryon stopping, and (2) the shorter lifetime of the fluid phase, which does not dilute the initial vorticity as much as it does at higher collision energies. At the same time, we did not observe a linear relation between the polarization and the ratio angular momentum/energy of the fireball.
A significant fraction of the produced originate from resonance decays. We have calculated the contribution to polarization stemming from the decay of polarized heavier states and particularly for the two leading contributing channels: 0 → γ and (1385) → π . It turns out that the contributions from 0 and (1385) alone change the resulting polarization by few percents, whereas extending the list of feed-down corrections up to mass 1670 MeV hyperons results in about 15% decrease of the resulting polarization.
It should be pointed out that we did not include posthadronization rescatterings in our calculation, which are likely to reduce the mean global polarization estimated in this work.
The calculated excitation function of the mean polarization reproduces the trend reported in the preliminary data from STAR [16,17], although the magnitude is about twice smaller, even without including the further reduction due to hadronic rescattering. As the uncertainty on the parameters of the model cannot compensate for this discrepancy, this result suggests that a revision of the initial state model, particularly of the initial longitudinal flow velocity profile is necessary (see e.g. [18]) as the present one provides an insufficient magnitude of shear longitudinal flow to approach the experimentally observed magnitude of the polarization.