Suppression of spin polarization as an indicator of QCD critical point

We study the impact of the QCD critical point (CP) on the spin polarization of Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}-hyperon generated by the thermal vorticity in viscous quark gluon plasma (QGP). The equations of the relativistic causal viscous hydrodynamics have been solved numerically in (3+1) dimensions to evaluate the thermal vorticity. The effects of the CP have been incorporated through the equation of state (EoS) and the scaling behavior of the transport coefficients. We observe a significant suppression of the spin polarization around mid-rapidity region induced by the CP and propose that such changed behaviour of the spin polarization can be used to detect the CP. To the best of our knowledge this is the first study on the effects of critical point on the spin polarization in QCD matter.


Introduction
Results from lattice quantum chromodynamics (QCD) and effective field theoretical models at non-zero temperature (T ) and baryon chemical potential (μ) reveal a rich and complex phase diagram [1]. While at high T and low μ (→ 0) the quark-hadron transition is a crossover, at low T and high μ the transition is believed to be first order in nature. Therefore, it is expected that between the crossover and the first order transition there exists a point in the μ − T plane called the Critical End Point or simply the Critical Point (CP) where the first order transition ends and crossover begins [2]. The location of CP is not yet known from first principles [3], however, phenomenological studies indicate its existence [4]. The investigations of the CP and the spin polarization in QCD matter are two frontier areas of research. In the present work a e-mail: sushant7557@gmail.com b e-mail: jane@vecc.gov.in (corresponding author) we introduce a connection between these two outstanding problems.
At present the general consensus is that the collisions of nuclei at the top Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) energies produce QGP with small μ and high T which reverts to hadronic phase via a crossover. The ongoing Beam Energy Scan-II program at RHIC and the upcoming Compressed Baryonic Matter experiment at Facility for Anti-proton and Ion Research (FAIR) and Nuclotron based Ion Collider fAcility (NICA) are planned to create system of quarks and gluons with different μ and T by varying the collision energies to explore the region close to the CP. The non-monotonic variation of the fluctuations in multiplicity with collision energy in the center of mass frame ( √ s N N ) [5], the change of sign of the fourth cumulant of order parameter with the variations of rapidity (rapidity scan) [6] and beam energy [7], the appearance of negative sign in the kurtosis of the order parameter fluctuation near the CP [8] are some of the proposed signals of the CP (see [9] for a review and references therein).
Relativistic hydrodynamics has been used to describe the space-time evolution of QGP and explain various experimental data quite successfully. One such crucial observable is the polarization of -hyperon [10] generated by the thermal vorticity and thermal shear during the hydrodynamic evolution of the QGP. Invigorating theoretical activities have been witnessed (see [11,12] for review) to understand various aspects of the polarization within the scope of hydrodynamics and transport models after the experimental measurement of the global polarization of the hyperon [10].
The magnetization of uncharged objects induced by mechanical rotation, called Barnett effect [13] and its inverse, that is, the rotation generated by varying magnetization, called the Einstein-de Haas effect [14] originate due to the conversion between spin (S) and orbital angular momentum (L) via spin-orbit coupling constrained by the conservation of total angular momentum J (= L + S). Similar coupling between L and S in QCD matter results in the spin polarization of particles. The initial orbital angular momentum (OAM) imparted by the spectators in non-central heavy ion collisions induce vortical motion locally and this induced vorticity polarize the quarks [15,16] which gets reflected in the spin polarization of final state hadrons. The coupling of fluid vorticity and quantum mechanical spin has been experimentally demonstrated for the first time in Ref. [17]. The vorticity in the fluid can be generated by other mechanisms (e.g. by the viscous stress as discussed in detail later). Hence, spin polarization of hadrons measured by STAR collaboration [10] has contributions from OAM and from other dynamical effects (e.g. viscosities, etc). The first contribution depends on the details of mechanism of transfer of initial OAM to the fluid and this is sensitive to the initial condition. The second contribution depends on the transport properties of the system and will be sensitive to the EoS. The goal of the present work is to understand the CP induced change in the vorticity and its consequences on the spin polarization of the hyperon. In other words if, CP μν ( μν ) is the vorticity in the presence (absence) of CP then what is the value of μν = ( CP μν − μν ) and the corresponding change on the spin polarization of -hyperon. We show in this work that as the CP is approached the local vorticity and hence the spin-polarization gets suppressed around the mid rapidity. Relativistic hydrodynamics can be applied to study the problem under consideration by choosing trajectories in the QCD phase diagram which are not too close to the CP but pass through the critical domain.
The presence of CP in the EoS affects the expansion of the system due to suppression of the sound wave and divergence of some of the transport coefficients [18]. This will affect the evolution of local vorticity and hence the -polarization through vorticity-spin coupling. Apart from polarization, the effect of CP on the separation of baryon and anti-baryon due to chiral vortical effect is another interesting facet [19].
Here we use natural unit c =h = k B = 1 where c is the speed of light in vacuum, h = 2πh is the Planck's constant and k B is the Boltzmann's constant. The signature metric for flat spacetime is taken as g μν = diag(1, −1, −1, −1).
The paper is organized as follows. In the next section we describe hydrodynamical model used in this work. In Sect. 3, we discuss the spin-polarization of -hyperon due to vorticity. Results are presented in Sects. 4 and 5 is devoted to summary and discussions.

Relativistic hydrodynamical model
We numerically solve (3+1)-dimensional relativistic viscous causal hydrodynamics using the algorithm detailed in Ref. [20]. The code contains the effect of CP through the EoS and the scaling behavior of the transport coefficients (see [21] for details). The initial condition and the EoS models used here to solve hydrodynamic equations have been extensively tested by reproducing the results available in Refs. [22] and [23] respectively. The CORNELIUS code [24] has been used to find the constant energy-density hyper surface. Our numerical results in the absence of CP have been contrasted with the known analytical results of Ref. [25] and with numerical results from other publicly available codes: AZHYDRO [26], MUSIC [27] and vHLLE [20]. The reliability of our code can be further appreciated by contrasting its output with the transverse momentum and rapidity dependence of various experimental observables, as we show in Sect. 4.
We have solved the following relativistic hydrodynamic equations numerically: where T μν is the energy-momentum tensor and N μ is the netbaryon number current. Here we work in the Landau frame of reference where the T μν and N μ are given by where is the bulk pressure, π μν is the shear-stress tensor which is symmetric, traceless and orthogonal to u μ , V μ (= 0, here) is the baryon diffusion 4-current and μν = g μν −u μ u ν . The viscous terms obey the following evolution equations, where · is defined as, N S and π μν N S are the Navier-Stokes limit of and π μν respectively, given by where θ = ∂ μ u μ . The coefficients of shear (η) and bulk (ζ ) viscosities are positive, i.e. η, ζ > 0. The hydrodynamical equations are solved in (τ, x, y, η s ) coordinates where, τ = √ t 2 − z 2 and η s = tanh −1 (z/t). The space-time evolution begins at time τ 0 . For lower energies the initial time, τ 0 is taken as the time required by the nuclei to pass through one other, ∼ 2R γ z v z (R is the radius of with v z is the velocity of the nucleus along z direction) and for higher energies ( √ s N N ≥ 62.4 GeV), τ 0 = 1 fm as shown in Table 1. The initial energy density profile at τ 0 is taken as: A symmetric rapidity profile, f (η s ), with the local energymomentum conservation puts a constraint on e(x, y) as shown in Ref. [22]. The energy deposited in the transverse plane, e(x, y) depends on the number of wounded nucleons per unit area which has been calculated by using the optical Glauber model for given impact parameter (b) at different √ s N N . The quantity, e(x, y) is related to the the number of wounded nucleons per unit area in the transverse plane, n A (x, y) and n B (x, y) of the colliding nuclei A and B respec- Here we consider Au+Au collisions at b = 5.6 fm for different √ s N N that has been found to correspond to 15-25% centrality [21]. The thickness function of the Au nucleus has been computed by assuming Woods-Saxon profile for nuclear density with nuclear radius, R = 6.37 fm, and surface thickness, δ = 0.535. The p+p inelastic cross-section, σ in N N ( √ s N N ), needed for the calculation of the number of wounded nucleons in the Glauber model has been taken from Refs. [28,29]. The initial velocity profile is taken as: The initial energy density and net baryon number density profiles have been computed with the parameters used in Ref. [22]. The viscous terms have been initialized with their corresponding Navier-Stokes limit. The EoS [23] employed here to solve the hydrodynamic equations reproduces the lattice QCD results at zero baryonic chemical potential. The parameters w, ρ and α 1 that appear in the linear mapping from Ising model to QCD in Ref. [23] have been fixed as w = 1, ρ = 2 and α 1 = 3.85 o . The other parameters are same as that of Ref. [23]. The transport coefficients are expected to diverge near the critical point following a scaling behavior [30]: where ξ(μ, T ) is the equilibrium correlation length, which is obtained through mapping QCD to 3D Ising model in the critical region. In the Ising model, ξ is computed by taking the derivative of equilibrium magnetization, M(r, h), with respect to the magnetic field, h, at fixed r = (T − T c )/T c as [30] where H 0 is a dimensionful parameter to get the correct dimensions of ξ . We shall take H 0 = 1 in our calculations and ∂ M ∂h r is obtained from the EoS model [23] using chain rule of differentiation. The extent of the critical domain in the μ − T plane is determined by the condition: ξ(μ, T ) = ξ 0 , where ξ 0 is taken as 1.75 fm. The possibility of divergent behavior is incorporated through the following expressions of the transport coefficients [30] Outside the critical region the values of the shear and bulk viscosities denoted by η 0 , ζ 0 respectively are chosen as [31,32]: The above parametrization is consistent with the estimates of the temperature-dependent specific shear and bulk viscosity extracted using Bayesian method [33] away from the critical region.
The dependence of the relaxation times appeared in Eqs. (4) and (5) on ξ are parameterized as: where τ 0 π and τ 0 are the relaxation times outside the critical region which are given by [31,32],

Thermal vorticity and spin polarization
The thermal vorticity at any space-time point of the fluid is given by [34][35][36]: where β μ = u μ /T . The -polarization is calculated by using the following expression of mean spin vector for a spin-1/2 particle with four-momentum p ν [37,38]: where m is the mass of the particle, μνρσ is the Levi-Civita tensor and n F is the Fermi-Dirac distribution. As we shall discuss later, the spin polarization is suppressed around mid rapidity in presence of CP. In Ref. [39], it has been shown that the contribution to spin polarization due to thermal shear around mid-rapidity is small. Hence, we have not considered the effect of thermal shear on spin polarization in this work for simplicity. Since the mass of the is much larger than the temperature range being considered in this study, we assume that 1 − n F ≈ 1 and n F ≈ n B , where n B is the Boltzmann distribution. Consequently the expression for the mean spin vector becomes In the rest frame of the particle, the spin vector is S * μ = (0, S * ), which is obtained by using the Lorentz transformation as: The mean spin averaged over the hypersurface is then given by [34], The net spin is obtained by integrating over azimuthal angle (0 ≤ φ < 2π ), rapidity (|y| < 1) and transverse momentum (0 < p T < 3 GeV) following the procedure of Ref. [40]. Finally the spin polarization of is given by,

Results
In this section, we present results through the following two subsections. First we test results from our newly developed computer code by contrasting with results from other codes available in the literature and experimental data. Then we discuss the effects of CP on spin polarization.

Testing the computer code
We present below a few test results from our hydrodynamic code. To test the code in (2+1)-dimensions, we compare the energy density obtained from the code with the Gubser's analytical solution [25] in (2+1) dimension in Fig. 1. We find that the results are in good agreement. There is no analytical result Fig. 1 The variation of ε with x. Results from the present work is compared with Gubser's analytical results in (2+1) dimension Fig. 2 The comparison of the rapidity distribution of π + from our code and the publicly available MUSIC code available in (3+1)-dimensions. Therefore, to test the code, we compare our result for the rapidity distribution of positively charged pion with the output from publicly available MUSIC code [41] without the resonance decays in Fig. 2.
Using the initial condition model, Eqs. (7) and (8) and EoS mentioned above, we further reproduce the charged particle pseudorapidity distribution for two colliding energies ( √ s N N = 19.6 GeV and 62.4 GeV) and different centralities in Fig. 3. To generate the plots of Fig. 3, we use the switching energy density ε sw = 0.3 GeV/fm 3 . We generate the constant energy density, ε sw , hypersurface henceforth denoted by ε , using the CORNELIUS code [24] which is then given as input to the UrQMD transport code [44] and generate 1000 events. It should be mentioned here that the width of the experimental distribution is slightly underestimated because we have used a single impact parameter and not an event-by-event simulation that would consist of a mixture of several impact parameters. We next reproduce the PHOBOS data on p T dependence of elliptic flow in 0-50% centrality of Au+Au collisions at √ s N N = 200 GeV [45] in Fig. 4. The variation of d N ch /dη with pseudorapidity with (EoS-CP) and without (noCP) the critical point in EoS are Fig. 3 The numerical results for d N ch /dη from our code are compared with experimental data of PHOBOS collaboration in different centralities [42,43]. Here EoS-CP and noCP denote results with and without the effects of CP respectively (see text) shown in Fig. 5 for colliding energies 14.5 GeV at impact parameter b = 4 fm (left), and 19.6 GeV at impact parameter b = 5.6 fm (right). The effect of the CP is found to be insignificant here.

Suppression of spin polarization by the CP
In this subsection, we will discuss the suppression of thermal vorticity and consequently of spin-polarization of -hyperon due to the presence of QCD CP. We will discuss the cases with zero as well as non-zero initial OAM.
We first show the critical region and the trajectories traced by the center of the fireball in the (μ − T ) plane at different √ s N N in Fig. 6. The black dot inside the critical domain (closed contour) indicates the location of the CP at (μ c , T c ) = (350MeV, 143.2MeV) [23] in this study. The trajectories have been calculated by solving the hydrodynamic equations with (denoted by EoS-CP) and without (denoted by noCP) the effects of CP. The trajectories for √ s N N =14.5 GeV and 19.6 GeV respectively passes through the critical domain. The trajectories for higher √ s N N remain outside the critical domain. We evaluate the thermal vorticity and subsequently the polarization of for system evolving along trajectories passing through both inside and outside the critical domain. As the transport code, UrQMD, does not include the spin effects, we stop the hydrodynamic evolution and evaluate the spin polarization on the constant energy density hypersurface ε = 0.3 GeV/fm 3 . The effect of the CP on the polarization is expected to be larger for √ s N N = 14.5 GeV as the trajectory for this case is closer to the CP compared to other values of √ s N N considered here.

Fig. 4
Our numerical result on p T dependence of elliptic flow, v 2 , compared to the experimental data from Ref. [45] Before presenting our results for non-zero initial OAM, some results with zero initial OAM are displayed first. The time evolution of the xη component of the thermal vorticity, averaged over the spatial coordinates and weighted by the energy density, with and without the effects of CP are shown in Fig. 7. Initially the system has zero vorticity. The vorticity generated by the viscous effects increase at first to attain some maximum value and then decreases subsequently. The evolution of the vorticity near the CP is affected by several competing factors like enhancement of various transport coefficients, slower expansion, changes in baroclinic torque, vortex stretching etc. Also, absorption of sound wave near the CP affects the expansion directly and the diverging nature of the transport coefficients reduces the vorticity as seen in Fig. 7. The observed suppression of the vorticity due It is intriguing to note that the CP not only reduces the polarization around mid-rapidity but also introduces strong qualitative changes in the slopes of the curves as shown in Fig. 9.
Finally, the global polarization is obtained as a function of √ s N N upon integration over p T , φ and y. The suppression of polarization is conspicuous for the trajectories passing through the critical domain at lower √ s N N (Fig. 10). The slope of the curve without CP is much steeper than the one with CP at lower √ s N N . The smaller value of polarization compared to experimental value [10] is because of the zero initial OAM. Inclusion of OAM through the non-zero initial value of the velocity profile [34] will enhance the magnitude of −P y , however, the difference in the rapidity distribution of polarization observed here with and without CP will still persist, as we discuss below. Now we introduce a non-zero OAM into the system through initial condition. The initial condition model that we use from Ref. [22] has been generalized to include a non-zero OAM in Ref. [46]. This is done by introducing a parameter, f , that takes value in the interval [0,1] and which controls the fraction of longitudinal momentum that can be attributed to the flow velocity. f = 0 corresponds to the Bjorken flow scenario. The assumption for the initial energy-momentum current in Ref. [46] can be achieved through the following choice of rest frame quantities: where p, ε, and u μ , respectively, denote the pressure, energy density, and fluid four flow-velocity. Also, y L (= f y CM ) denotes the local longitudinal rapidity variable and y CM is the rapidity in local center-of-mass frame. The components of the initial energy-momentum tensor then take the following form: which is consistent with the assumption of Ref. [46]. The trace of the energy-momentum tensor is given by T μ μ = T τ τ − τ 2 T ηη = 0. Hence, our choice for the rest frame quantities, Eq. (14), does not violate the condition for the trace of the energy- GeV in 20-60% centrality [47]. Vertical lines are the statistical uncertainties only 20-60% centrality [47]. Vertical lines are the statistical uncertainties only momentum tensor. It was stated in Ref. [46] that the parameter f has negligible effects on most of the global observables such as the pseudorapidity distributions, particle yields, and elliptic flow. We have checked this and conclude the same. By setting the value of f = 0.2, impact parameter, b = 8.7 fm, and using the same set of values for other parameters of the IC model for Au+Au collision at √ s N N = 200 GeV, described earlier, we compute the negative y-component of global polarization, −P y , of -hyperon for 20-60% centrality. We get −P y =0.254% where the corresponding experimental value from STAR is 0.277± 0.040 (stat). By using the updated PDG value of α , we get −P y = 0.243± 0.035% (stat) [47]. Having fixed the parameter f with the help of global polarization data, we set out to see if we can reproduce the measurements of local polarization. We show the pseudorapidity and transverse momentum dependence in the 20-60% centrality in Figs. 11 and 12 respectively. It is to be stressed that we have used EoS with CP for the results displayed in Figs. 11 and 12 for √ s NN = 200 GeV. The After qualitatively understanding the CP induced effects, we make prediction for the spin polarization as the critical point is approached. For this we choose √ s N N = 14.5 GeV which is close to the CP in our model study. We first tune the parameter f so as to reproduce the STAR measurement of global polarization at this energy [10]. The values of f = 0.45 with CP and f = 0.53 without CP corresponds to the same value of −P y =0.92% of -hyperon for b = 5.6 which reproduces the experimental value measured by the STAR collaboration for 20-60% centrality. After tuning f we evaluate the rapidity distribution of P y . We plot the rapidity distribution of −P y in the left panel of Fig. 13 for these values of f . A suppression of about 25% in −P y is observed at mid-rapidity. The changes in other observables like p T spectra, elliptic flow, d N ch /dy is at most 8% on the surface ε [21]. We also compute the derivative of P y with respect to rapidity and plot as a function of rapidity in the right panel of Fig. 13. We observe a slight positive slope for d P y /dy at mid-rapidity opposite to the case when there is no critical point. We cannot confirm the sign of the slope by further approaching the critical point due to the limitations of the EoS model which is valid for μ B < 450 MeV. The observed suppression of spin polarization at √ s N N = 14.5 GeV has resulted even though the fireball is 100 MeV away from the critical point along the μ B axis. We expect the effect to get enhanced on further approach toward the critical point.
The sensitivity of the spin polarization on the EoS can be understood from the following expression for the spin polarization in the rest frame of -hyperon at any point on ε [11]: where ω is the vorticity (= ∇ × v), γ is the Lorentz factor and A denotes the acceleration of the fluid element. In a nutshell, spin polarization depend on the gradients of temperature and curl of flow-velocity (= ω). These gradients in turn depend on the expansion dynamics of the system and the expansion is strongly influenced by the speed of sound through EoS.
Since the sound mode gets suppressed at the critical point, the system undergoes a slow expansion consequently the gradients of temperature and flow-velocity develop slowly. This should then result into a suppression of spin-polarization as confirmed by our simulations.

Summary and discussions
The thermal vorticity and consequently the polarization of the hyperon for different colliding energies have been estimated and found to be suppressed as the CP is approached. It is well-known that the local vorticity of the fluid couples with the quantum mechanical spin of the particles and polarize them. We have evaluated the spin polarization of -hyperon with and without the effects of CP and found a strong change in the spin polarization around mid-rapidity as the system approaches the CP. There are various physical processes which collectively contribute to the suppression. Although we have solved the fully relativistic equation to estimate the vorticity, we consider below the evolution equation for vorticity ( ω = ∇ × v) for a compressible fluid with constant ζ and η in the non-relativistic limit because the contributions from various terms appear more clearly in this limit, Here ρ denotes the density of fluid and θ = ∇ · v. As θ is a measure of the expansion of the system, a larger expansion results in smaller vorticity as suggested by the negative sign of the term θ ω, the vortex stretching term, ω · ∇ v enhances it; the baroclinic torque described by 1 ρ 2 ∇ρ × ∇ p generates vorticity by creating shear due to different acceleration confronted by the fluid with different densities. The term proportional to ∇ 2 ω is responsible for spatial diffusion of vorticity, the diffusion coefficient being η ρ . The term of particular interest is proportional to ∇ρ × ∇θ which suggests that the vorticity dissipates if there is a gradient in expansion rate for fluid cells i.e. the fluid cells having less density and expanding faster will oppose the vorticity of the denser fluid cells expanding slowly. The strength of this effect is proportional to ζ + 1 3 η . It is clear that the suppression of vorticity and hence polarization, is a combined effect of the absorption of sound wave and the enhancement of various transport coefficients in presence of CP. The drastic qualitative and quantitative changes induced by CP in the rapidity distribution of P y can be used to detect the CP experimentally. It is important to mention at this point that the effects of CP on the p T spectra of the hadrons and on the p T and y dependence of directed and elliptic flow are found to be small [21].
Some comments on the application of hydrodynamics near the CP are in order here. Near the CP, the fluctuating modes do not relax faster than the timescale of changes in slow/conserved variables due to which the local thermal equilibrium is not maintained making hydrodynamics inapplicable. The validity of the hydrodynamics can, however, be extended by adding a scalar variable representing the slow non-hydrodynamic modes connected to the relaxation rate of the critical fluctuation (see [48] and [49] for details). It has been explicitly shown that the modes associated with the scalar variable lags behind the hydrodynamic modes resulting in back reactions on the hydrodynamic variables [50]. Further, it has been demonstrated in Ref. [50] that the back reaction has negligible effects on the hydrodynamic variables. In view of this, the results presented in this work will be useful in detecting the CP. The addition of slow nonhydrodynamic modes related to the relaxation of the critical fluctuation will be considered in future studies. Moreover, we may also recall that if a system is not too close to CP then hydrodynamics can still be applied in a domain around the CP [51].