A study of vorticity formation in high energy nuclear collisions

We present a quantitative study of vorticity formation in peripheral ultrarelativistic heavy-ion collisions at sNN=200\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sqrt{s}_\mathrm{NN}}= 200$$\end{document} GeV by using the ECHO-QGP numerical code, implementing relativistic dissipative hydrodynamics in the causal Israel–Stewart framework in 3+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3+1$$\end{document} dimensions with an initial Bjorken flow profile. We consider different definitions of vorticity which are relevant in relativistic hydrodynamics. After demonstrating the excellent capabilities of our code, which proves to be able to reproduce Gubser flow up to 8 fm/c, we show that, with the initial conditions needed to reproduce the measured directed flow in peripheral collisions corresponding to an average impact parameter b=11.6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=11.6$$\end{document} fm and with the Bjorken flow profile for a viscous Quark Gluon Plasma with η/s=0.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta /s=0.1$$\end{document} fixed, a vorticity of the order of some 10-2c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-2}\,c$$\end{document}/fm can develop at freeze-out. The ensuing 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} baryons does not exceed 1.4 % at midrapidity. We show that the amount of developed directed flow is sensitive to both the initial angular momentum of the plasma and its viscosity.


Introduction
The hydrodynamical model has by now become a paradigm for the study of the QCD plasma formed in nuclear collisions at ultrarelativistic energies. There has been a considerable advance in hydrodynamics modeling and calculations of these collisions over the last decade. Numerical simulations in 2 + 1 D [1] and in 3 + 1 D [2][3][4][5][6][7] including viscous corrections are becoming the new standard in this field, and existing codes are also able to handle initial state fluctuations. a e-mail: becattini@fi.infn.it An interesting issue is the possible formation of vorticity in peripheral collisions [8][9][10]. Indeed, the presence of vorticity may provide information as regards the (mean) initial state of the hydrodynamical evolution which cannot be achieved otherwise, and it is related to the onset of peculiar physics in the plasma at high temperature, such as the chiral vortical effect [11]. Furthermore, it has been shown that vorticity gives rise to polarization of particles in the final state, so that e.g. baryon polarization -if measurable -can be used to detect it [12,13]. Finally, as we will show, a numerical calculation of vorticity can be used to make stringent tests of numerical codes, as the T -vorticity (see Sect. 2 for the definition) is expected to vanish throughout under special initial conditions in the ideal case.
Lately, vorticity has been the subject of investigations in Refs. [9,10] with peculiar initial conditions in cartesian coordinates, ideal fluid approximation and isochronous freezeout. Instead, in this work, we calculate different kinds of vorticity with our 3 + 1 D ECHO-QGP 1 code [3], including dissipative relativistic hydrodynamics in the Israel-Stewart formulation with Bjorken initial conditions for the flow (i.e. with u x = u y = u η = 0), henceforth denoted BIC. It should be pointed out from the very beginning that the purpose of this work is to make a general assessment of vorticity at top RHIC energy and not to provide a precision fit to all the available data. Therefore, our calculations do not take into account effects such as viscous corrections to the particle distribution at the freeze-out and initial state fluctuations, that is, we use smooth initial conditions obtained averaging over many events.
We will use the relativistic notation with repeated indices assumed to be summed over, however, contractions of indices will be sometimes denoted with dots, e.g. u·T ·u ≡ u μ T μν u ν . The covariant derivative is denoted d μ (hence d λ g μν = 0), the exterior derivative by d, whereas ∂ μ is the ordinary derivative.

Vorticities in relativistic hydrodynamics
Unlike in classical hydrodynamics, where vorticity is the curl of the velocity field v, several vorticities can be defined in relativistic hydrodynamics which can be useful in different applications (see also the review [14]).

The kinematical vorticity
The kinematical vorticity is defined as where u is the four-velocity field. This tensor includes both the acceleration A and the relativistic extension of the angular velocity pseudo-vector ω μ in the usual decomposition of an antisymmetric tensor field into a polar and pseudo-vector fields: where μνρσ is the Levi-Civita symbol. Using the transverse (to u) projector: and with the usual definition of the orthogonal derivative where D = u α d α , it is convenient to define also a transverse kinematical vorticity as Using the above definition in the decomposition (2) it can be shown that that is, ω is the tensor formed with the angular velocity vector only. As we will show in the next subsection, only ω shares the "conservation" property of the classical vorticity for an ideal barotropic fluid.

The T -vorticity
This is defined as and it is particularly useful for a relativistic uncharged fluid, such as the QCD plasma formed in nuclear collisions at very high energy. This is because from the basic thermodynamic relations when the temperature is the only independent thermodynamic variable, the ideal relativistic equation of motion (ε + p)A μ = ∇ μ p can be recast in the simple form (see e.g. [15]): Equation (6) is also known as the Carter-Lichnerowicz equation [14] for an ideal uncharged fluid and it entails conservation properties which do not hold for the kinematical vorticity. This can be better seen in the language of differential forms, rewriting the definition of the T -vorticity as the exterior derivative of a vector field (1-form) T u, that is = d(T u). Indeed, Eq. (6) implies -through the Cartan identity -that the Lie derivative of along the vector field u vanishes, that is, because is itself the external derivative of the vector field T u and dd = 0. Equation (7) states that the T -vorticity is conserved along the flow and, thus, if it vanishes at an initial time it will remain so at all times. This can be made more apparent by expanding the Lie derivative definition in components: The above equation is in fact a differential equation for precisely showing that if = 0 at the initial time then ≡ 0. Thereby, the T -vorticity has the same property as the classical vorticity for an ideal barotropic fluid, such as the Kelvin circulation theorem, so the integral of over a surface enclosed by a circuit comoving with the fluid will be a constant.
One can write the relation between T -vorticity and kinematical vorticity by expanding the definition (5): implying that the double-transverse projection of , μρ νσ Hence, the tensor ω shares the same conservation properties of , namely it vanishes at all times if it is vanishing at the initial time. Conversely, the mixed projection of the kinematical vorticity, does not. It then follows that for an ideal uncharged fluid with ω = 0 at the initial time, the kinematical vorticity is simply

The thermal vorticity
This is defined as [13]: where β is the temperature four-vector. This vector is defined as (1/T )u once a four-velocity u, that is a hydrodynamical frame, is introduced, but it can also be taken as a primordial quantity to define a velocity through u ≡ β/ β 2 [16]. The thermal vorticity features two important properties: it is adimensional in natural units (in cartesian coordinates) and it is the actual constant vorticity at the global equilibrium with rotation [17] for a relativistic system, where β is a Killing vector field whose expression in Minkowski space-time is β μ = b μ + μν x ν , b and being constant. In this case the magnitude of thermal vorticity is -with the natural constants restored -simplyhω/k B T where ω is a constant angular velocity. In general (replacing ω with the classical vorticity defined as the curl of a proper velocity field) it can be readily realized that the adimensional thermal vorticity is a tiny number for most hydrodynamical systems, though it can be significant for the plasma formed in relativistic nuclear collisions. Furthermore, the thermal vorticity is responsible for the local polarization of particles in the fluid according to the formula [12] μ (x, p) = − which applies to spin 1/2 fermions, n F being the Fermi-Dirac-Juttner distribution function, Similarly to the previous subsection, one can readily obtain the relation between T -vorticity and thermal vorticity: Again, the double-transverse projection of is proportional to the one of : whereas the mixed projection turns out to be, using Eq. (13), Again, for an ideal uncharged fluid with ω = 0 at the initial time, by using the equations of motion (6), one sees that the above projection is just A ν /T and that the thermal vorticity is simply A common feature of the kinematical and thermal vorticity is that their purely spatial components can be non-vanishing if the acceleration and velocity field are non-parallel, even though velocity is vanishing at the beginning.

High energy nuclear collisions
In nuclear collisions at very large energy, the QCD plasma is an almost uncharged fluid. Therefore, according to previous section's arguments, in the ideal fluid approximation, if the transversely projected vorticity tensor ω initially vanishes, so will the transverse projection and and the kinematical and thermal vorticities will be given by the formulas (9) and (14), respectively. Indeed, the T -vorticity will vanish throughout because also its longitudinal projectionvanishes according to Eq. (6). This is precisely what happens for the usually assumed BIC for the flow at τ 0 , that is u x = u y = u η = 0, where one has ω = 0 at the beginning as it can be readily realized from the definition (1). On the other hand, for a viscous uncharged fluid, transverse vorticities can develop even if they are zero at the beginning.
It should be noted, though, that even if the space-space components (x, y, η indices) of the kinematical vorticity tensor vanish at the initial Bjorken time τ 0 , they can develop at later times even for an ideal fluid if the spatial parts of the In the full longitudinally boost invariant Bjorken picture, that is, u η = 0 throughout the fluid evolution, in the ideal case, as ω = 0, the only allowed components of the kinematical vorticity are ω τ x , ω τ y and ω xy from the first Eq. (2). The ω xy component, at η = 0, because of the reflection symmetry (see Fig. 1) in both the x and the y axes, can be different from zero but it ought to change sign by moving clockwise (or counterclockwise) to the neighboring quadrant of the x y plane; for central collisions it simply vanishes.
However, in the viscous case, more components of the vorticities can be non-vanishing. Furthermore, in more realistic 3 + 1 D hydrodynamical calculations, a non-vanishing u η can develop because of the asymmetries of the initial energy density in the x − η and y − η planes at finite impact parameter. The asymmetry is essential to reproduce the observed directed flow coefficient v 1 (y) in a 3+1 D ideal hydrodynamic calculation with BIC, as shown by Bozek and Wyskiel [18], and gives the plasma a total angular momentum, as will be discussed later on.
In this work, we calculate the vorticities, and especially the thermal vorticity by using basically the same parametrization of the initial conditions in Ref. [18]. Those initial conditions are a modification of the usual BIC to take into account that the plasma, in peripheral collisions, has a relatively large angular momentum (see Appendix A). They are a minimal modifications of the BIC in that the initial flow velocity Bjorken components are still zero, but the energy density longitudinal profile is changed and no longer symmetric by the reflection η → −η. They are summarized hereinafter. Given the usual thickness function expression: where n 0 = 0.1693 fm −3 , δ = 0.535 fm, and R = 6.38 fm are the nuclear density, the width and the radius of the nuclear Fermi distribution respectively, the following functions are defined: where σ is the inelastic NN cross section, A the mass number of the colliding nuclei, and where x T = (x, y) is the vector of the transverse plane coordinates and b is the impact parameter vector, connecting the centers of the two nuclei. In our conventional cartesian reference frame, the b vector is oriented along the positive x axis and the two nuclei have initial momentum along the z axis (whence the reaction plane is the xz plane) and their momenta are directed so as to make the initial total angular momentum oriented along the negative y axis (see Fig. 1). The wounded nucleons weight function W N is then defined: where Finally, the initial proper energy density distribution is assumed to be where the total weight function W (x, y, η) is defined as and In Eq. (21) n BC (x, y) is the mean number of binary collisions: and α is the collision hardness parameter, which can vary between 0 and 1. This parametrization, and especially the chosen forms of the functions f ± , are certainly not unique as a given angular momentum can be imparted to the plasma in infinitely many ways. Nevertheless, as has been mentioned, it proved to reproduce correctly the directed flow in a 3+1 D hydrodynamical calculation of peripheral Au-Au collisions at high energy [18], thus we took it as a good starting point. A variation of this initial condition will be briefly discussed in Sect. 7. Besides, the parametrization (20) essentially respects the causality constraint that the plasma cannot extend beyond η = y beam . Indeed, at √ s NN = 200 GeV y beam 5.36 while the 3σ point in the gaussian profile in Eq. (22) lies at η = η flat /2 + 3σ η 4.4.
The free parameters have been chosen following Ref. [19], where they were adjusted to reproduce the data in Au-Au collisions at √ s NN = 200 GeV. They are reported in Table 1. We have run the ECHO-QGP code in both the ideal and the viscous modes with the parameters reported in Table 1 [20]. The impact parameter value b = 11.57 was chosen as, in the optical Glauber model, it corresponds to the mean value of the 40-80 % centrality class (9.49 < b < 13.42 fm [21]) used by the STAR experiment for the directed flow measurement in Ref. [22]. The initial flow velocities u x , u y , u η were set to zero, according to BIC. The freeze-out hypersurfaceisothermal at T fo = 130 MeV -is determined with the methods described in Refs. [3,23].

Qualification of the ECHO-QGP code
To show that our code is well suited to model the evolution of the matter produced in heavy-ion collisions and hence to carry out our study on the development of vorticity in such an environment, we have performed two calculations, referring to an ideal and viscous scenario, respectively, providing a very stringent numerical test. Before describing these tests, it should be pointed out that the vorticities components are to be calculated in Bjorken coordinates, whose metric tensor is g μν = diag(1, −1, −1, −τ 2 ), hence they do not all have the same dimension nor are they adimensional, as is desirable (except the thermal vorticity, as has been emphasized in Sect. 2). For a proper comparison it is better to use the orthonormal basis, which involves a factor τ when the η components are considered. Moreover, the cumulative contribution of all components is well described by the invariant modulus, which, for a generic antisymmetric tensor A μν , is Furthermore, we have always rescaled the T -vorticity by 1/T 2 in order to have an adimensional number. Since the T -vorticity has always been determined at the isothermal freeze-out, in order to get its actual magnitude, one just needs to multiply it by T 2 fo .

T -vorticity for an ideal fluid
Since the fluid is assumed to be uncharged and the initial Tvorticity is vanishing with the BIC, it should be vanishing throughout, according to the discussion in Sect. 2. However, the discretization of the hydrodynamical equations entails a numerical error, thus the smallness of in an ideal run is a gauge of the quality of the computing method. In Fig. 2 we show the mean of the absolute values of the six independent Bjorken components at the freeze-out hypersurface, of the Tvorticity divided by T 2 to make it adimensional, as a function of the grid resolution (the boundaries in x, y, η being fixed). 2 As is expected, the normalized T -vorticity decreases as the resolution improves. Because of Eq. (13), the residual value at our best spatial resolution of 0.15 fm can be taken as a numerical error for later calculations of the thermal vorticity (Fig. 3).

Gubser flow
A very useful test for the validation of a numerical code of relativistic dissipative hydrodynamics is the explicit solution of Israel-Stewart theory of a Bjorken flow with an azimuthally symmetric radial expansion [24][25][26][27], the so-called Gubser flow. Indeed, this solution provides a highly non-trivial theoretical benchmark. 2 It should be pointed out that, throughout this work, by mean values of the vorticities we mean simple averages of the Bjorken components (possibly rescaled by 1/τ ) over the freeze-out hypersurface without geometrical cell weighting. Therefore, the plotted mean values have no physical meaning and they should be taken as descriptive numbers which are related to the global features of vorticity components at the freeze-out.
For the sake of clarity, we briefly summarize the main steps leading to the analytical solution, to be compared with the numerical computation. In the case of a conformal fluid, with p = /3 EOS, the invariance for scale transformations sets the terms entering the second-order viscous hydrodynamic equations. The additional requirements of azimuthal and longitudinal-boost invariance constrain the solution of the hydrodynamic equations, which has to be invariant under SO(3) q ⊗ SO(1, 1) ⊗ Z 2 transformations. To start with, one defines a modified space-time metric as follows (with the usual Bjorken coordinates, η being the space-time rapidity): which can be viewed as a rescaling of the metric tensor: It can be shown that dŝ 2 is the invariant space-time interval of d S 3 ⊗ R, where d S 3 is the three-dimensional de Sitter space and R refers to the rapidity coordinate. It is then convenient to perform a coordinate transformation (q is an arbitrary parameter setting an energy scale for the solution once one goes back to physical dimensionful coordinates) after which the rescaled space-time element dŝ 2 reads The full symmetry of the problem is now manifest. SO(1, 1) and Z 2 refer to the usual invariance for longitudinal boosts and η → −η inversion, while SO(3) q reflects the spherical symmetry of the rescaled metric tensor in the new coordinates. In Gubser coordinates the fluid is at rest: The corresponding flow in Minkowski space can be obtained taking into account both the rescaling of the metric and the change of coordinates wherex μ = (ρ, θ, φ, η) and x μ = (τ, r, φ, η). Other quantities such as the temperature or the viscous tensors require the solution of the following set of hydrodynamic equations (their most general form actually admits further terms that were derived for a system of massless particles in Refs. [28,29]), valid for the case of a conformal fluid with ε = 3 p: In the case of the Gubser flow in Eq. (28), due to the traceless and transverse conditionsπ μ μ = 0 andû μπ μ ν = 0, one has simply to solve the two equations (π ηη ≡π ηη /ŝT ) and (η/ŝ = η/s, the ratio being dimensionless) The solution can then be mapped back to Minkowski space through the formulas In Fig. 4 we show the comparison between the Gubser analytical solution and our numerical computation for the temperature T and the components π x x , π xy , and π ηη of the viscous stress tensor, respectively, at different times. The initial energy density profile was taken from the exact Gubser solution at the time τ = 1 fm/c. The simulation has been performed with a grid of 0.025 fm in space and 0.001 fm in time. The shear viscosity to entropy density ratio was set to η/s = 0.2, while the shear relaxation time is τ R = 5η/(ε + p). The energy scale is set to q = 1 fm −1 . As it can be seen, the agreement is excellent up to late times.

T -vorticity for a viscous fluid
Unlike for an ideal uncharged fluid, T -vorticity can be generated in a viscous uncharged fluid even if it is initially vanishing. Thus, the T -vorticity can be used as a tool to estimate the numerical viscosity of the code in the ideal mode by extrapolating the viscous runs.
A comment is in order here. In general, in addition to standard truncation errors due to finite-difference interpolations, all shock-capturing upwind schemes are known to introduce numerical approximations that behave roughly as a dissipative effects, especially in the simplified solution to the Riemann problems at cell interfaces [30]. It is therefore important to check whether the code is not introducing, for a given resolution, numerical errors which are larger than the effects induced by the physics. We refer to the global numerical errors generically as numerical viscosity.
We have thus calculated the T -vorticity for different physical viscosities (in fact η/s ratios), in order to provide an upper bound for the numerical viscosity of ECHO-QGP in the ideal mode. The mean value of the T -vorticity is shown in Fig. 5 and its extrapolation to zero occurs when |η/s| 0.002 which is a very satisfactory value, comparable with the one obtained in Ref. [4]. The good performance is due to the use of high-order reconstruction methods that are able to compensate for the highly diffusive two-wave Riemann solver employed [3].

Directed flow, angular momentum, and thermal vorticity
With the initial conditions reported at the end of Sect. 3 we have calculated the directed flow of pions (both charged states) at the freeze-out and compared it with the STAR data for charged particles collected in the centrality interval 40-80 % [22]. Directed flow is an important observable for several reasons. Recently, it has been studied at lower energy [31] with a hybrid fluid-transport model (see also Ref. [32]). At √ s NN = 200 GeV, it has been calculated with an ideal 3 + 1 D hydro code first by Bozek and Wyskiel [18]. Herein, we extend the calculation to the viscous regime.
The amount of generated directed flow at the freeze-out depends, of course, on the initial conditions, particularly on the parameter η m (see Sect. 3), as shown in Fig. 6. The directed flow also depends on η/s as shown in Fig. 7 and could then be used to measure the viscosity of the QCD plasma along with other azimuthal anisotropy coefficients. It should be pointed out that, apparently, the directed flow can be reproduced by our hydrodynamical calculation only for −3 < y < 3.
The dependence of v 1 (y) on η m and η/s makes it possible to adjust the η m parameter for a given η/s value. This adjustment cannot be properly called a precision fit because, as we have mentioned in the Introduction, several effects in the  [22] comparison between data and calculations have been deliberately neglected in this work. However, since our aim was to obtain a somewhat realistic evaluation of the vorticities, we have chosen the value of η m for which we obtain the best agreement between our calculated pion v 1 (y) and the measured for charged particles in the central rapidity region. For the fixed value η/s = 0.1 (approximately twice the conjectured universal lower bound) the corresponding best value of η m turns out to be 2.0 (see Fig. 8).
It is worth discussing more in detail an interesting relationship between the value of the parameter η m and that of a conserved physical quantity, the angular momentum of the plasma, which, for BIC is given by the integral (see Appendix A for the derivation): Since η m controls the asymmetry of the energy density distribution in the η − x plane, one expects that J y will vary as a function of η m . Indeed, if the energy density profile is symmetric in η, the integral in Eq. (32) vanishes. Yet, for any finite η m = 0, the profile (20) is not symmetric and J y = 0 (looking at the definition of f + and f − it can be realized that only in the limit η m → ∞ the energy density profile becomes symmetric). The dependence of the angular momentum on η m with all the initial parameters kept fixed is shown in Fig. 9.
For the value η m = 2.0 it turns out to be around 3.18 × 10 3 inh units. It is also interesting to estimate an upper bound on the angular momentum of the plasma by evaluating the angular momentum of the overlap region of the two colliding nuclei. This can be done by trying to extend the simple formula for two sharp spheres. In our conventional reference frame, the initial angular momentum of the nuclear overlap region is directed along the y axis with negative value and can be written as where T ± are the thickness functions like in Eq. (18) and w(x, y) = min(n(x + b/2, y, 0), n(x − b/2, y, 0)) max(n(x + b/2, y, 0), n(x − b/2, y, 0)) is the function which extends the simple product of two θ functions used for the overlap of two sharp spheres. Note that theω(x, y) is 1 for full overlap (b = 0) and implies a vanishing angular momentum for very large b (see Fig. 10) (see also Ref. [33]). At b = 11.57 fm the above angular momentum is about 3.58 × 10 3 inh units. This means that, with the current parametrization of the initial conditions, for that impact parameter about 89 % of the angular momentum is retained by the hydrodynamical plasma, while the rest is possibly taken away by the corona particles.
With the final set of parameters, we have calculated the thermal vorticity . As has been mentioned in Sect. 2, this vorticity is adimensional in cartesian coordinates) and it is constant at global thermodynamical equilibrium [17], e.g. It can be seen that the generated amount of thermal vorticity has some non-trivial dependence on the viscosity. Particularly, as it is apparent from Fig. 12, the xη componentwhich is directed along the initial angular momentum -has a non-vanishing mean value whose magnitude significantly increases with increasing viscosity. Its map at the freeze-out, for a fixed value of the y coordinate y = 0, is shown in Fig. 13 where it can be seen that it attains a top (negative) value of about 0.05 corresponding to a kinematical vorticity, at the freeze-out temperature of 130 MeV, of about 0.033 c/fm 10 22 s −1 . In this respect, the Quark Gluon Plasma would be the fluid with the highest vorticity ever made in a terrestrial laboratory. However, the mean value of this component at the same value of η/s = 0.1 is of the order of 5.4 × 10 −3 , that is about ten times less than its peak value,  Fig. 12. This mean thermal vorticity is consistently lower than the one estimated in Ref. [13] (about 0.05) with the model described in Refs. [9,10], implying an initial non-vanishing transverse kinematical and thermal vorticity . This reflects in a quite low value of the polarization of baryons, as will be shown in the next section.

Polarization
As has been mentioned in the Introduction, vorticity can result in the polarization of particles in the final state. The relation between the polarization vector of a spin 1/2 particle and thermal vorticity in a relativistic fluid was derived in Ref. [12] and reads where n F is the Fermi-Dirac-distribution function (12) and the integration is over the freeze-out hypersurface . The interesting feature of this relation is that it makes it possible to obtain an indirect measurement of the mean thermal vorticity at the freeze-out by measuring the polarization of some hadron. For instance, the polarization of baryons, as is well known, can be determined with the analysis of the angular distribution of its decay products, because of parity violation. The polarization pattern depends on the momentum of the decaying particle, as it is apparent from Eq. (34). Equation (34) makes sense only if the components of the integrand are Minkowskian, as an integrated vector field yields a vector only if the tangent spaces are the same at each point. Before summing over the freeze-out hypersurface we have then transformed the components of the thermal vorticity from Bjorken coordinates to Minkowskian by using the known rules. The polarization vector ( p) thus obtained is the one in the collision frame. However, the polarization vector which is measurable is the one in the decaying particle In Fig. 14 we show the polarization vector components, as well as its modulus, as a function of the transverse momentum p T for p z = 0 expected under the assumptions of local thermodynamical equilibrium for the spin degrees of freedom maintained till kinetic freeze-out. It can be seen that the polarization vector has quite an assorted pattern, with an overall magnitude (see Fig. 14, panel a) hardly exceeding 1 % at momenta around 4 GeV. As expected, the y component is predominantly negative, oriented along the initial angular momentum vector and a magnitude of the order of 0.1 %. Indeed, the main contribution to the polarization stems from the longitudinal component z 0 , with a maximum and minimum along the bisector | p x | = |p y |.
The obtained polarization values are -as expected -consistently smaller than those estimated in Ref. [13] (of the order of several percent with a top value of 8-9 %) with the already mentioned initial conditions used in Refs. [9,10]. This is a consequence of the much lower value of the implied thermal vorticity, as discussed in the previous section. Also, the y 0 pattern is remarkably different, with different location of maxima and minima.

Conclusions, discussion and outlook
To summarize, we have calculated the vorticities developed in peripheral (b = 11.6 fm) nuclear collisions at √ s NN = 200 GeV (b = 11.6 fm) with the most commonly used initial conditions in the Bjorken hydrodynamical scheme, by using the code ECHO-QGP implementing second-order, causal, relativistic dissipative hydrodynamics. An extensive testing of the high accuracy and very low numerical diffusion properties of the code has been carried out, followed by longtime simulations (up to τ = 8 fm/c) of the so-called viscous Gubser flow, a stringent test of numerical implementations of Israel-Stewart theory in Bjorken coordinates.
We have found that the magnitude of the 1/τ x − η component of the thermal vorticity at freeze-out can be as large as 5 × 10 −2 , and yet its mean value is not large enough to produce a polarization of hyperons much larger than (36) compared with STAR data [22] both the shear viscosity and the longitudinal energy density profile asymmetry parameter η m which in turn governs the amount of initial angular momentum retained by the plasma. The fact that in 3 + 1 D the plasma needs to have an initial angular momentum in order to reproduce the observed directed flow raises the question whether the Bjorken initial condition u η = 0 is a compelling one or, instead, the same angular momentum can be obtained with a non-trivial u η and with a suitable change of the energy density profile. For testing purposes, we have run ECHO-QGP with the initial profile u η = 1 τ tanh Ax sinh(y beam − |η|), which meets the causality constraint (see Appendix B). It is found that the directed flow is very sensitive to an initial u η . For a small positive value of the parameter A = 5 × 10 −4 fm −1 corresponding to a J y = 3.32 × 10 3 , keeping all other parameters fixed, the directed flow exhibits two slight wiggles around midrapidity (see Fig. 15) which are not seen in the data. For a very small negative value of the parameter A = −5×10 −4 fm −1 , corresponding to J y = 3.08×10 3 , the directed flow increases while approximately keeping the same shape as for A = 0 around midrapidity. However, more detailed studies are needed to determine whether a non-vanishing initial flow velocity is compatible with the experimental observables. We plan to extend this kind of calculation to different centralities, different energies and with initial state fluctuations in order to determine the possibly best conditions for vorticity formation in relativistic nuclear collisions.