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

TL;DR

This paper analytically investigates baryon polarization in heavy ion collisions using exact hydrodynamical solutions, highlighting the role of rotation and vorticity in generating polarization during the quark-gluon plasma evolution.

Contribution

It provides the first analytical calculations of polarization four-vectors in heavy ion collisions using exact hydrodynamical solutions, including a rotating and accelerating expansion model.

Findings

Polarization depends on the fluid's vorticity and rotation.

Analytic solutions show the necessity of rotation for non-zero polarization.

Results illustrate how different hydrodynamical flows influence baryon polarization.

Abstract

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.

Equations

$$f(x,p)\propto \frac{1}{\exp \left(\frac{p_{\mu }{u}^{\mu }(x)}{T(x)}-\frac{\mu (x)}{T(x)}\right)+1},$$

$$\displaystyle\langle S(x,p)\rangle^{\mu}=\frac{1}{8m}\big{(}1{-}f(x,p)\big{)}\varepsilon^{\mu\nu\rho\sigma}p_{\sigma}\partial_{\nu}\beta_{\rho},$$

$$E\frac{\mathrm{d}N}{{\mathrm{d}}^3\mathbf{p}}=\int {\mathrm{d}}^3{\Sigma }_{\mu }(x)p^{\mu }f(x,p).$$

$$⟨S(p){⟩}^{\mu }=\frac{\int {\mathrm{d}}^3{\Sigma }_{\nu }{p}^{\nu }f(x,p)⟨S(x,p){⟩}^{\mu }}{\int {\mathrm{d}}^3{\Sigma }_{\nu }{p}^{\nu }f(x,p)},$$

$$\displaystyle\int\mathrm{d}^{d}\mathbf{r}\,f(\mathbf{r})g(\mathbf{r})\approx f(\mathbf{R}_{0})g(\mathbf{R}_{0})\sqrt{\frac{(2\pi)^{d}}{\det\mathbf{M}}},\qquad\begin{array}[]{ll}\textnormal{where}&\quad\mathbf{M}_{ij}=\partial_{i}\partial_{j}g(\mathbf{r})\big{|}_{{\mathbf{r}}={\mathbf{R}_{0}}},\\ \textnormal{and}&\quad\partial_{k}g(\mathbf{R}_{0})=0,\end{array}$$

$$\displaystyle f(x,p)\ll 1\quad\quad\Rightarrow\quad\quad f(x^{\mu},p^{\mu})=\frac{g}{(2\pi\hbar)^{d}}\exp\big{(}\frac{\mu(x)}{T(x)}{-}\frac{p_{\mu}u^{\mu}}{T(x)}\big{)}.$$

$$\displaystyle f(x^{\mu},p^{\mu})=C_{0}\exp\big{(}{-}p_{\mu}\beta^{\mu}(x)\big{)},\qquad\textnormal{where}\quad\beta^{\mu}(x)=\frac{u^{\mu}(x)}{T(x)},\quad\textnormal{and}\quad C_{0}=\frac{g}{(2\pi\hbar)^{d}}.$$

$$⟨S(x,p){⟩}^{\mu }=\frac{1}{8m}{\epsilon }^{\mu \nu \rho \sigma }{p}_{\sigma }{\partial }_{\nu }{\beta }_{\rho },$$

$$\displaystyle\langle S(p)\rangle^{\mu}\approx\frac{1}{8m}\varepsilon^{\mu\nu\rho\sigma}p_{\sigma}\partial_{\nu}\beta_{\rho}\Big{|}_{\mathbf{r}{=}\mathbf{R}_{0}},$$

$$(\epsilon +p)u^{\nu }{\partial }_{\nu }{u}^{\mu }$$

$$(\epsilon +p){\partial }_{\mu }{u}^{\mu }$$

$$n{\partial }_{\mu }{u}^{\mu }$$

$$u^{\mu }$$

$$n=n_0{\left(\frac{{\tau }_0}{\tau }\right)}^d\mathcal{V}(S),T=T_0{\left(\frac{{\tau }_0}{\tau }\right)}^{d/\kappa }\frac{1}{\mathcal{V}(S)},$$

$$S\equiv S(S_x,S_y,S_z),\text{where}{S}_x\equiv \frac{r_x^2}{{\dot{X}}_0^2{t}^2},S_y\equiv \frac{r_y^2}{{\dot{X}}_0^2{t}^2},S_z\equiv \frac{r_z^2}{{\dot{X}}_0^2{t}^2},\text{for example:}S=\frac{r_x^2}{{\dot{X}}_0^2{t}^2}+\frac{r_y^2}{{\dot{Y}}_0^2{t}^2}+\frac{r_z^2}{{\dot{Z}}_0^2{t}^2}.$$

$${\mathrm{d}}^3{\Sigma }_{\mu }=\frac{1}{t(\mathbf{r})}\left(\begin{array}{c}t(\mathbf{r})\\ \mathbf{r}\end{array}\right){\mathrm{d}}^3\mathbf{r}\Rightarrow E\frac{\mathrm{d}N}{{\mathrm{d}}^3\mathbf{p}}=C_0\int {\mathrm{d}}^3\mathbf{r}\frac{Et(\mathbf{r})-\mathbf{pr}}{t(\mathbf{r})}\exp \left(-\frac{Et(\mathbf{r})-\mathbf{p}\mathbf{r}}{T_0}\right).$$

$$\displaystyle\partial_{k}\frac{Et{-}\mathbf{p}{\mathbf{r}}}{T_{0}}\Big{|}_{\mathbf{r}{=}\mathbf{R}_{0}}\stackrel{{\scriptstyle!}}{{=}}0\quad\Rightarrow\quad\mathbf{R}_{0}=\frac{\tau_{0}}{m}\mathbf{p}.\qquad\qquad M_{kl}\equiv-\partial_{k}\partial_{l}\frac{Et-\mathbf{p}{\mathbf{r}}}{T_{0}}\Big{|}_{\mathbf{r}=\mathbf{R}_{0}}=\frac{m}{T_{0}\tau_{0}}\left({\delta_{kl}-\frac{p_{k}p_{l}}{E^{2}}}\right).$$

$$\det \mathbf{M}=\frac{m^2}{E^2}{\left(\frac{m}{T_0{\tau }_0}\right)}^3\Rightarrow E\frac{\mathrm{d}N}{{\mathrm{d}}^3\mathbf{p}}=\frac{n_0}{4}\sqrt{\frac{\pi {\tau }_0^3}{m{T}_0^3}}\exp \left(-\frac{{\tau }_{0}m}{T_0}\right).$$

$${\partial }_{\nu }{\beta }_{\rho }={\partial }_{\nu }\left(\frac{r_{\rho }}{\sqrt{{\tau }_0^2+r^2}{T}_0}\right)=\frac{g_{\nu \rho }}{\sqrt{{\tau }_0^2+r^2}{T}_0}+\frac{r_{\nu }{r}_{\rho }}{({\tau }_0^2+r^2{)}^{3/2}{T}_0}.$$

$$\displaystyle\langle S(p)\rangle^{0}=\frac{1}{8mT_{0}}\varepsilon^{0ikl}p_{l}\partial_{i}\beta_{k}\bigg{|}_{\mathbf{r}=\mathbf{R}_{0}}=\frac{1}{8mT_{0}}\varepsilon_{ikl}p_{l}\left({\frac{g_{ik}}{\sqrt{\tau_{0}^{2}{+}r^{2}}T_{0}}+\frac{r_{i}r_{k}}{(\tau_{0}^{2}{+}r^{2})^{3/2}T_{0}}}\right)\Bigg{|}_{\mathbf{r}{=}\mathbf{R}_{0}}=0,$$

$$\displaystyle\langle S(p)\rangle^{i}=\frac{1}{8mT_{0}}\bigg{(}-\varepsilon_{ikl}p_{l}\partial_{0}\beta_{k}+\varepsilon_{ikl}p_{l}\partial_{k}\beta_{0}-\varepsilon_{ikl}p_{0}\partial_{k}\beta_{l}\bigg{)}\Bigg{|}_{\mathbf{r}=\mathbf{R}_{0}}=0.$$

$$⟨S(p){⟩}^{\mu }=\left(\begin{array}{c}0\\ 0\end{array}\right),$$

$$\mathbf{v}$$

$$\frac{u^{\mu }}{T}\equiv {\beta }^{\mu }=a^{\mu }+F^{\mu \nu }{x}_{\nu }+(x^{\nu }{b}_{\nu })x^{\mu }-\frac{x^{\nu }{x}_{\nu }}{2}{b}^{\mu },$$

$$a^{\mu }=\frac{{\rho }_0^2}{2T_0{\tau }_0^2}\left(\begin{array}{c}1\\ 0\end{array}\right),b^{\mu }=\frac{1}{T_0{\tau }_0^2}\left(\begin{array}{c}1\\ 0\end{array}\right),F_{0k}=F_{k0}=F_{00}=0,F_{kl}={\epsilon }_{klm}\frac{{\Omega }_m}{2T_0}.$$

$$\displaystyle E\frac{\mathrm{d}N}{\mathrm{d}^{3}\mathbf{p}}=C_{0}\int\mathrm{d}^{3}\mathbf{r}\left({E{-}\frac{\mathbf{p}\mathbf{r}}{\sqrt{\tau_{0}^{2}{+}r^{2}}}}\right)\exp\left\{{-\frac{E(2r^{2}{+}\tau_{0}^{2}{+}\rho_{0}^{2}){-}2\sqrt{\tau_{0}^{2}{+}r^{2}}\mathbf{p}\mathbf{r}{-}\tau_{0}^{2}\mathbf{r}(\mathbf{p}{\times}\mbox{\boldmath${\Omega}$})}{T_{0}\tau_{0}^{2}}}\right\}.$$

$$\displaystyle\textnormal{for $\mathbf{R}_{0}$}:\qquad\nabla\Big{\{}{-}\frac{1}{T_{0}\tau_{0}^{2}}\Big{(}E(2r^{2}{+}\tau_{0}^{2}{+}\rho_{0}^{2})-2\sqrt{\tau_{0}^{2}{+}r^{2}}\mathbf{r}\mathbf{p}-\tau_{0}^{2}\mathbf{r}(\mathbf{p}{\times}\mbox{\boldmath${\Omega}$})\Big{)}\Big{\}}\Big{|}_{\mathbf{r}=\mathbf{R}_{0}}\stackrel{{\scriptstyle!}}{{=}}0,$$

$$\displaystyle M_{kl}=\partial_{k}\partial_{l}\Big{\{}\frac{1}{T_{0}\tau_{0}^{2}}\Big{(}E(2r^{2}{+}\tau_{0}^{2}{+}\rho_{0}^{2})-2\sqrt{\tau_{0}^{2}{+}r^{2}}\mathbf{r}\mathbf{p}-\tau_{0}^{2}\mathbf{r}(\mathbf{p}{\times}\mbox{\boldmath${\Omega}$})\Big{)}\Big{\}}\Big{|}_{\mathbf{r}=\mathbf{R}_{0}}.$$

$$\displaystyle\mathbf{R}_{0}=\frac{\tau_{0}}{2p}\sqrt{\frac{E{-}m}{2m}}\sqrt{\tau_{0}^{2}(\hat{\mathbf{p}}{\times}\mbox{\boldmath${\Omega}$})^{2}(E{-}m)^{2}+4p^{2}}\cdot\hat{\mathbf{p}}+\tau_{0}^{2}\frac{E{-}m}{2p}\cdot\hat{\mathbf{p}}{\times}\mbox{\boldmath${\Omega}$}.$$

$$\det M_{kl}=\frac{32m^2}{T_0^3{\tau }_0^6}(E+m)p.$$

$$\displaystyle E\frac{\mathrm{d}N}{\mathrm{d}^{3}\mathbf{p}}\propto\sqrt{\frac{\pi^{3}T_{0}^{3}\tau_{0}^{3}}{32p(m{+}E)}}\exp\left({-\frac{E_{\textnormal{eff}}}{T_{0}}}\right),\quad\textnormal{with}\quad E_{\textnormal{eff}}=m{+}\frac{\rho_{0}^{2}E}{\tau_{0}^{2}}{+}\frac{\tau_{0}^{2}}{4}(\mbox{\boldmath${\Omega}$}^{2}{-}(\hat{\mathbf{p}}\mbox{\boldmath${\Omega}$})^{2})\left({E{-}m}\right).$$

$$\displaystyle E\frac{\mathrm{d}N}{\mathrm{d}^{3}\mathbf{p}}\propto\sqrt{\frac{\pi^{3}T_{0}^{3}\tau_{0}^{3}}{32p(m{+}E)}}\exp\left({-\frac{E}{T_{\textnormal{eff}}}}\right),\quad\textnormal{with}\quad T_{\textnormal{eff}}=\frac{T_{0}}{\frac{m}{E}{+}\frac{\rho_{0}^{2}}{\tau_{0}^{2}}{+}\frac{\tau_{0}^{2}}{4}(\mbox{\boldmath${\Omega}$}^{2}{-}(\hat{\mathbf{p}}\mbox{\boldmath${\Omega}$})^{2})\left({1{-}\frac{m}{E}}\right)}.$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date \TitlePolarized baryon production in heavy ion collisions: an analytic hydrodynamical study \AuthorB. Boldizsár 1\orcidA, M. I. Nagy 1 and M. Csanád 1 \AuthorNamesB. Boldizsár, M. I. Nagy, M. Csanád

\corresCorrespondence: nmarci@elte.hu

1 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 $\Lambda$ baryons at the STAR experiment STAR:2017ckg ; Adam:2018ivw 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. Becattini:2013fla .

Although many numerical hydrodynamical models do indeed predict non-zero polarization of produced spin 1/2 particles Csernai:2014nva ; Xie:2016fjj ; Karpenko:2016jyx ; Xie:2017upb , 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:

[TABLE]

where $p_{\mu}$ is the four-momentum of the produced particle, and $u^{\mu}(x)$, $\mu(x)$ and $T(x)$ are the four-velocity, the chemical potential, and the temperature field of the fluid, respectively.

Assuming local thermal equlibrium for the spin degrees of freedom, for the space-time– and momentum-dependent polarization four-vector $\langle S(x,p)\rangle^{\mu}$ of the produced particles the following formula is given in Ref. Becattini:2013fla :

[TABLE]

where $m$ is the mass of the investigated particle, and the inverse temperature field $\beta^{\mu}{=}u^{\mu}/T(x)$ is introduced. Here $\varepsilon^{\mu\nu\rho\sigma}$ is the totally antisymmetric Levi-Civita-symbol; the $\varepsilon^{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 Csorgo:2003rt ; Csorgo:2003ry and a rotating and accelerating solution (first reported in Ref. Nagy:2009eq , then in a different context in Hatta:2014gqa ). 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.

2 Basic equations and assumptions

We use the $c{=}1$ notation. Let us denote the space-time coordinate by $x^{\mu}{\equiv}(t,\mathbf{r})$, and the Minkowskian metric tensor by $g^{\mu\nu}{=}\mathrm{diag}(1,{-}1,{-}1,{-}1)$. The convention for the Levi-Civita symbol is $\varepsilon^{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^{\mu}}_{\mu}=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^{\mu}{=}\gamma(1,\mathbf{v})$, where $\gamma{=}\sqrt{1{-}v^{2}}$ is the Lorentz factor. The velocity three-vector is then $\mathbf{v}{=}u^{k}/u^{0}$. With $p^{\mu}$ we denote the four-momentum of a produced particle; we also use the three-momentum $\mathbf{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=\sqrt{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

[TABLE]

Here $\mathrm{d}^{3}\Sigma_{\nu}$ is the 3-dimensional vectorial integration measure of the freeze-out hypersurface; the appearance of which is the so-called Cooper-Frye prescription Cooper:1974mv 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. Becattini:2013fla for the polarization of spin 1/2 particles, as written up in Eq. (2), may be utilized for any given $\beta^{\mu}{=}u^{\mu}/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 $\langle S(x,p)\rangle^{\mu}$ field over the freeze-out hypersurface. The formula to be analyzed further, that is, that for the observed polarization $\langle S(p)\rangle^{\mu}$ of particles with momentum $p$, thus becomes

[TABLE]

as written up e.g. in Karpenko:2016jyx . 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(\mathbf{r})g(\mathbf{r})$, where $f(\mathbf{r})$ is a slowly changing function, while $g(\mathbf{r})$ has a unique and sharp maximum; then the integral can be calculated with a Gaussian approximation as

[TABLE]

that is, $\mathbf{R}_{0}$ is the location of the unique maximum of $g(\mathbf{r})$ and $\mathbf{M}$ is the second derivative matrix.

Another assumption concerns the expression of $\langle S(x,p)\rangle^{\mu}$, 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:

[TABLE]

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 $\mu/T$ factor can (and usually is) neglected; we use this approximation here111 The vanishing of $\mu$ 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 $\mu/T=\mathrm{const}$ instead of $\mu/T=0$; if $\mu\neq 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 $\int\mathrm{d}p\,f(x,p)=n(x)$ is evident. For example, if one chooses an ultra-relativistic ideal gas, with $p=nT$, $\varepsilon=\kappa p$, with $\kappa=d$ as EoS, one has $\frac{g}{(2\pi\hbar)^{d}}e^{\mu/T}=\frac{n}{4\pi T^{3}}$. Indeed, in the solutions discussed below, $\mu/T=$const is satisfied, which means $n\propto T^{d}$, which is the well-known condition for an adiabatic expansion. . With this we have

[TABLE]

If the Maxwell-Boltzmann approximation is justified, it means that $f(x,p)\ll 1$ indeed, and then also Eqs. (2) and (4) become simpler:

[TABLE]

and in the saddle-point approximation, the polarization of particles with momentum $p$ becomes simply

[TABLE]

since in the saddle-point approximation, in the numerator of Eq. (4), $\langle S(x,p)\rangle^{\mu}$ can be considered the ,,smooth” function, and the determinant factors cancel.

3 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

[TABLE]

and we specify the simple $\varepsilon=\kappa p$ equation of state here. (The notations: $\varepsilon$, $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 $\mu=0$). So the expressions for $n$ that we recapitulate for the solutions can be regarded as supplemental to the solutions that work for $\mu=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 Florkowski:2018fap , along with some numerical calculations of how this modified hydrodynamical picture affects final state polarization Florkowski:2019qdp . We do not investigate this possibility here; we restrict ourselves to the simple and well-known basic equations witten up above.

3.1 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. Csorgo:2003rt . This solution has the following velocity, particle density and temperature fields:

[TABLE]

where $\tau=\sqrt{t^{2}-\mathbf{r}^{2}}$, and $\kappa$ is the inverse square speed of sound (constant in the case of this exact solution). The $\kappa{=}3$ case corresponds to ultrarelativistic ideal gas, $\kappa{=}3/2$ corresponds to a non-relativistic gas; however, this solution is valid for any arbitrary constant $\kappa$ value222We note that a more general class of solutions is possible Csorgo:2003rt ; Csorgo:2003ry ; Csanad:2014dpa in which the temperature and density fields are supplemented with an arbitrary $\mathcal{V}$ function of a ,,scaling variable” $S$:

$n=n_{0}\left(\frac{\tau_{0}}{\tau}\right)^{d}\mathcal{V}(S),\qquad T=T_{0}\left(\frac{\tau_{0}}{\tau}\right)^{d/\kappa}\frac{1}{\mathcal{V}(S)},$

and the $S$ variable is any function of $S_{x}$, $S_{y}$, and $S_{z}$:

$S\equiv S(S_{x},S_{y},S_{z}),\quad\textnormal{where}\quad S_{x}\equiv\frac{r_{x}^{2}}{\dot{X}_{0}^{2}t^{2}},\quad S_{y}\equiv\frac{r_{y}^{2}}{\dot{X}_{0}^{2}t^{2}},\quad S_{z}\equiv\frac{r_{z}^{2}}{\dot{X}_{0}^{2}t^{2}},\quad\textnormal{for example:}\quad S=\frac{r_{x}^{2}}{\dot{X}_{0}^{2}t^{2}}+\frac{r_{y}^{2}}{\dot{Y}_{0}^{2}t^{2}}+\frac{r_{z}^{2}}{\dot{Z}_{0}^{2}t^{2}}.$

Here $\dot{X}_{0}$, $\dot{Y}_{0}$ and $\dot{Z}_{0}$ are arbitrary constants. In the given example, the $S=\mathrm{const}$ surfaces are ellipsoids, and $\dot{X}_{0}$, $\dot{Y}_{0}$, $\dot{Z}_{0}$ are time derivatives of the principal axes of them.

To calculate the polarization four-vector, as of now we investigate the simplest case, the spherical symmetric expansion. For the freeze-out hypersurface the $\tau{=}\tau_{0}{=}\mathrm{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 $\mathbf{r}$ coordinate three-vector, and the time coordinate on the hypersurface is $t(\mathbf{r}){\equiv}\sqrt{\tau_{0}^{2}{+}\mathbf{r}^{2}}$. The integration measure and the resulting expression for the Cooper–Frye formula can then be written as

[TABLE]

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 ($\mathbf{R}_{0}$) as well as the second derivative matrix $M_{kl}$ is calculated as

[TABLE]

With this we can get an approximation for the invariant single-particle momentum distribution:

[TABLE]

The formula is independent of momentum. This was expected because this hydrodynamical solution (in the $\mathcal{V}(S){=}1$ case) is boost invariant.

To use (11) to determine the polarization four-vector in the hydrodynamical solution of the Hubble-flow, first we give the expression for the $\partial_{\nu}\beta_{\rho}$ derivative:

[TABLE]

Then for the time component we get:

[TABLE]

as $\varepsilon^{0ikl}$ is antisymmetric whereas $g_{ik}$ and $r_{i}r_{k}$ are symmetric to the change in the $i\leftrightarrow k$ indices.

Similarly for the spatial coordinates:

[TABLE]

In conclusion, the polarization four-vector in the spherical symmetric Hubble-flow is

[TABLE]

which is consistent with our expectations.

3.2 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. Nagy:2009eq . This solution has the following velocity, temperature and particle density profiles:

[TABLE]

where $\rho_{0}$ and $\tau_{0}$ are arbitrary parameters and ${\Omega}$ is an arbitrary angular velocity three-vector that indicates the axis and magnitude of rotation. The $\rho_{0}$ parameter tells about the initial spatial extent of the expanding matter, however, the $\tau_{0}$ parameter is just there for the sake of consistency of physical units; in this way, the unit of ${\Omega}$ is $c$/fm, as it should be for an angular velocity-like quantity333Here we changed the notation of Ref. Nagy:2009eq . The rather unfortunate $\mathbf{B}$ notation used there is now written as \tau_{0}^{2}\mbox{\boldmath{\Omega}}., and $T_{0}$ is a temperature constant. In the case of \mbox{\boldmath{\Omega}}{=}0, we recover an acceleratingly expanding but non-rotating spherically symmetric solution.

It is convenient to write up this solution with the following notation:

[TABLE]

To calculate final state observables, we choose the constant proper time ($\tau_{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 $\mathbf{r}=0$ center of the expanding matter). We use the notation introduced in Eq. (13) 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:

[TABLE]

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 ($\mathbf{R}_{0}$) and the second derivative matrix at the saddle point:

[TABLE]

We leave the detailed calculations to Appendix .1; the results are the following. The $\mathbf{R}_{0}$ saddle point (for a given $\mathbf{p}$ momentum) is in the plane spanned by the $\mathbf{p}$ and \mathbf{p}{\times}\mbox{\boldmath{\Omega}} vectors. In the following we use the $\hat{\mathbf{p}}\equiv\mathbf{p}/p$ notation for the unit vector pointing in the direction of $\mathbf{p}$. For the saddle point we get

[TABLE]

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

[TABLE]

Using this result, we get the invariant single-particle momentum distribution444This has not yet been calculated for this hydrodynamical solution. as

[TABLE]

Equivalently, by defining a “local slope” $T_{\textnormal{eff}}$, the result can be expressed as

[TABLE]

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. (21), then substitute it into the expression of the polarization, Eq. (11). The result is

[TABLE]

(The second term was cancelled owing to the symmetry of $g_{\nu\rho}$ and the antisymmetry of $\varepsilon^{\mu\nu\rho\sigma}$, and $x^{\mu}$ is understood as the four-coordinate of the freeze-out hypersurface whose three-coordinate is the $\mathbf{r}=\mathbf{R}_{0}$ three-vector). Remembering the expression of the introduced $F^{\mu\nu}$ tensor and $b^{\mu}$ vector from Eq. (21), in particular that $F^{0k}=0$, and $b^{k}=0$, we get the following expressions for the the time-like and space-like components:

[TABLE]

Summarizing this result, the polarization four-vector for the investigated rotating and accelerating expanding solution is the following:

[TABLE]

In the case of \mbox{\boldmath{\Omega}}=0, there is no rotation, and we get $\langle S(p)\rangle^{\mu}{=}0$. In this model thus polarization is very transparently connected to the presence of rotation.

It is useful to transform the polarization four-vector into the rest frame of the particle. The result is555 The Lorentz matrix performing this boost transformation is the following (in usual 1+3 dimensional block matrix notation):

${\Lambda^{\mu}}_{\nu}=\begin{pmatrix}\cosh\chi&-\hat{p}_{l}\sinh\chi\\ -\hat{p}_{k}\sinh\chi&\delta_{kl}{+}(\cosh\chi{-}1)\hat{p}_{k}\hat{p}_{l}\end{pmatrix}=\frac{1}{m}\begin{pmatrix}E&-p_{l}\\ -p_{k}&m\delta_{kl}+\frac{E{-}m}{p^{2}}p_{k}p_{l}\end{pmatrix},$

where $E$ and $p$ could be parametrized with the velocity parameter $\chi$ as $E=m\cosh\chi$ and $p=m\sinh\chi$, respectively. It indeed can be checked that this matrix takes the $(E,\mathbf{p})$ four-momentum vector into $(m,\mathbf{0})$, as it should., with (r.f. standing for “rest frame”):

[TABLE]

We can also compute the helicity of the produced spin 1/2 particles in this solution from this formula (the $\mathbf{S}$ polarization vector is taken in the laboratory frame):

[TABLE]

4 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. Karpenko:2016jyx ). 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 $\Lambda$ baryon ($m_{\Lambda}c^{2}=1115$ MeV) was chosen. For the sake of this illustration, we chose a moderate value for the magnitude of the ${\Omega}$ vector as |\mbox{\boldmath{\Omega}}|=0.1\,c/\mathrm{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 $\mathbf{p}$; see Eq. (34). 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 $\mathbf{p}\mathbf{S}$ 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.

5 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 gradient have on the observable final state polarization of baryons produced in heavy-ion collisions. Such studies have the potential of a better understanding of what phenomenological implications can polarization measurements (such as recently done by the STAR experiment STAR:2017ckg ) have on the properties (such as the equation of state) of the strongly coupled Quark Gluon Plasma produced in heavy ion collisions.

Acknowledgements

\funding

Our research has been partially supported by the Hungarian NKIFH grants No. FK-123842 and FK-123959, the Hungarian EFOP 3.6.1-16-2016-00001 project. M. Csanád and M. Nagy was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences and the ÚNKP New National Excellence Program of the Hungarian Ministry of Human Capacities.

.1 Additional calculations

Here we discuss some additional calculations used in Section 3.2 pertaining to the case of rotating and accelerating solution.

For a given momentum $\mathbf{p}$, the position of the saddle point $\mathbf{R}_{0}$ (to be applied in the approximate calculation of the momentum spectrum and the polarization) was written up in Eq. (26); we provide some additional details of the derivation of that formula here. The defining equation was Eq. (24), of which the following equation for $\mathbf{R}_{0}$ is obtained:

[TABLE]

where $R_{0}^{2}{\equiv}\mathbf{R}_{0}\mathbf{R}_{0}$. From this equation one readily sees that $\mathbf{R}_{0}$ must be a linear combination of $\mathbf{p}$ and the \mathbf{p}{\times}\mbox{\boldmath{\Omega}} vector. We substitute this assumption into the equation above. We note that $\mathbf{p}$ and \mathbf{p}{\times}\mbox{\boldmath{\Omega}} are orthogonal to each other, which leads to some intermediate simplifications, as well as enables us to rearrange the obtained condition into the following form:

[TABLE]

where we temporarily introduced the A{\equiv}\sqrt{\tau_{0}^{2}{+}\alpha^{2}p^{2}+\beta^{2}\tau_{0}^{4}(p^{2}\Omega^{2}{-}(\mathbf{p}\mbox{\boldmath{\Omega}})^{2})} notation. Because of the orthogonality of $\mathbf{p}$ and \mathbf{p}{\times}\mbox{\boldmath{\Omega}}, both sides here have to vanish identically, from which we get

[TABLE]

One divides these equations to obtain a simple relation, the substituting back one gets a quadratic equation for $\beta$, the solution of which is

[TABLE]

where we used the $E^{2}{=}p^{2}{+}m^{2}$ relation. To find $\alpha$ we substitute this back into the expression of $A$:

[TABLE]

Using the above expression of $\beta$ (with the yet undetermined sign) we get $1{-}4p^{2}\beta^{2}=-\frac{2m}{p^{2}}(m{\pm}E)$, and see that the expression for $\alpha$ will be valid only in the case when $1{-}4\beta^{2}p^{2}>0$, thus conclude that the bottom sign is the proper choice. We thus arrive at the following expressions:

[TABLE]

From these formulas the expression of $\mathbf{R}_{0}$ shown in Eq. (26) readily follows. The other ingredient in the saddle-point integration necessary for getting the momentum spectrum is the determinant of the second derivative matrix of the source function. Here we outline the main steps of the derivation of Eq. (27). From Eq. (25) the second derivative matrix itself turns out to be

[TABLE]

where we use the notation $A$ as above. We should use the expression of $\mathbf{R}_{0}$ as calculated above.

The determinant of this $\mathbf{M}$ matrix is the product of its eigenvalues. In our case the particular spatial directions are: $\mathbf{p}$, \mathbf{p}{\times}\mbox{\boldmath{\Omega}}, and the vector orthogonal to both these, that is, \mathbf{p}{\times}(\mathbf{p}{\times}\mbox{\boldmath{\Omega}}). One recognizes that the vector \mathbf{p}{\times}(\mathbf{p}{\times}\mbox{\boldmath{\Omega}}) is an eigenvector of the $\mathbf{M}$ second derivative matrix:

[TABLE]

The corresponding eigenvalue is thus $1/\beta$. Owing to the symmetric nature of $\mathbf{M}$, the other two eigenvectors must be in the orthogonal complementer subspace of this vector, so they are linear combinations of $\mathbf{p}$ and \mathbf{p}{\times}\mbox{\boldmath{\Omega}}. Let us thus look for these eigenvectors in the form $\mathbf{a}{=}\mu\mathbf{p}{+}\nu\mathbf{r}$, with yet to be determined $\mu$ and $\nu$ coefficients. Substituting this expression, we get

[TABLE]

where $\lambda$ is the eigenvalue (the values of which we are looking for). By substituting the assumed form of $\mathbf{a}$ and inferring the components of this equation in the $\mathbf{p}$ and \mathbf{p}{\times}\mbox{\boldmath{\Omega}} directions, we get the following equation for the $\mu$ and $\nu$ coefficients:

[TABLE]

We immediately infer the product of the two $\lambda_{1,2}$ eigenvalues as the determinant of this $2{\times}2$ matrix. Taking the third eigenvalue (calculated above) into account, after some simplifications, we indeed get the following expression for the determinant of the $\mathbf{M}$ matrix (the expression we used in Eq. (27)):

[TABLE]