(3+1)-dimensional dissipative relativistic fluid dynamics at non-zero net baryon density

TL;DR

This paper introduces BEShydro, a (3+1)-dimensional relativistic hydrodynamic code that models high-energy heavy-ion collisions with finite baryon density, including dissipative effects, and provides extensive testing protocols.

Contribution

The paper presents a new modular (3+1)-D hydrodynamic code for simulating baryon-rich QCD matter with detailed dissipative effects and comprehensive validation tests.

Findings

BEShydro accurately models baryon diffusion and viscous effects.

The code's modularity allows selective inclusion of physical effects.

Validation tests demonstrate high precision in baryon transport simulations.

Abstract

Heavy-ion collisions at center-of-mass energies between 1 and 100 GeV/nucleon are essential to understand the phase diagram of QCD and search for its critical point. At these energies the net baryon density of the system can be high, and simulating its evolution becomes an indispensable part of theoretical modeling. We here present the (3+1)-dimensional diffusive relativistic hydrodynamic code BEShydro which solves the equations of motion of second-order Denicol-Niemi-Molnar-Rischke (DNMR) theory, including bulk and shear viscous currents and baryon diffusion currents. BEShydro features a modular structure that allows to easily turn on and off baryon evolution and different dissipative effects and thus to study their physical effects on the dynamical evolution individually. An extensive set of test protocols for the code, including several novel tests of the precision of baryon transport that can also be used to test other such codes, is documented here and supplied as a permanent part of the code package.

Equations

$$d_{\mu }{T}^{\mu \nu }$$

$$d_{\mu }{N}^{\mu }$$

$$\Gamma^{\mu}_{\alpha\beta}\equiv\frac{1}{2}g^{\mu\nu}\bigl{(}\partial_{\beta}g_{\alpha\nu}+\partial_{\alpha}g_{\nu\beta}-\partial_{\nu}g_{\alpha\beta}\bigr{)}=\Gamma^{\mu}_{\beta\alpha}\;.$$

$$T^{\mu \nu }$$

$$N^{\mu }$$

$$T^{\mu \nu }{u}_{\nu }=\mathcal{E}{u}^{\mu },$$

$$\mathcal{E}=u_{\mu }{T}^{\mu \nu }{u}_{\nu },\mathcal{N}=u_{\mu }{N}^{\mu }.$$

$$D\mathcal{N}$$

$$D\mathcal{E}$$

$$(\mathcal{E}+{\mathcal{P}}_0+\Pi )D{u}^{\mu }$$

$${\tau }_{\Pi }D\Pi +\Pi$$

$${\tau }_n(Dn{)}^{⟨\mu ⟩}+n^{\mu }$$

$${\tau }_{\pi }(D\pi {)}^{⟨\mu \nu ⟩}+{\pi }^{\mu \nu }$$

$${\Pi }_{\mathrm{NS}}$$

$$n_{\mathrm{NS}}^{\mu }$$

$${\pi }_{\mathrm{NS}}^{\mu \nu }$$

$${\tau }_{\Pi }D\Pi +\Pi$$

$${\tau }_{n}D{n}^{\mu }+n^{\mu }$$

$${\tau }_{\pi }D{\pi }^{\mu \nu }+{\pi }^{\mu \nu }$$

$$(Dn{)}^{⟨\mu ⟩}={\Delta }_{\nu }^{\mu }D{n}^{\nu }=D{n}^{\mu }+u^{\mu }{n}^{\nu }D{u}_{\nu },$$

$$\tau =\sqrt{t^2-z^2},{\eta }_s=\frac{1}{2}\ln \left(\frac{t+z}{t-z}\right).$$

$$g^{\mu\nu}=\mathrm{diag}\bigl{(}1,-1,-1,-1/\tau^{2}\bigr{)}\;,$$

$${\Gamma }_{\tau \eta }^{\eta }={\Gamma }_{\eta \tau }^{\eta }=\frac{1}{\tau },{\Gamma }_{\eta \eta }^{\tau }=\tau .$$

$${\partial }_{\mu }{T}^{\mu \tau }$$

$${\partial }_{\mu }{T}^{\mu x}$$

$${\partial }_{\mu }{T}^{\mu y}$$

$${\partial }_{\mu }{T}^{\mu \eta }$$

$${\partial }_{\mu }{N}^{\mu }$$

$$d\Pi$$

$$d{n}^{\mu }$$

$$d{\pi }^{\mu \nu }$$

$$I_{\Pi }$$

$$I_n^{\mu }$$

$$I_{\pi }^{\mu \nu }$$

$$I_1^{\mu }=u^{\mu }{n}^{\nu }D{u}_{\nu },I_2^{\mu }=n^{\mu }\theta ,I_3^{\mu }=n_{\nu }{\omega }^{\nu \mu },I_4^{\mu }=n_{\nu }{\sigma }^{\nu \mu };$$

$$I_1^{\mu \nu }=\left(u^{\mu }{\pi }^{\nu \lambda }+u^{\nu }{\pi }^{\mu \lambda }\right)D{u}_{\lambda },I_2^{\mu \nu }=\theta {\pi }^{\mu \nu },I_3^{\mu \nu }={\omega }_{\lambda }^{\mu }{\pi }^{\lambda \nu }+{\omega }_{\lambda }^{\nu }{\pi }^{\lambda \mu },$$

$$I_4^{\mu \nu }=\frac{1}{2}\left({\pi }^{\mu \lambda }{\sigma }_{\lambda }^{\nu }+{\pi }^{\nu \lambda }{\sigma }_{\lambda }^{\mu }\right)-\frac{1}{3}{\Delta }^{\mu \nu }{\pi }^{\alpha \beta }{\sigma }_{\beta \alpha }.$$

$$\frac{\eta }{{\tau }_{\pi }}$$

$$\frac{{\delta }_{\pi \pi }}{{\tau }_{\pi }}=\frac{4}{3},\frac{{\tau }_{\pi \pi }}{{\tau }_{\pi }}$$

$$\bar{\eta }=\frac{\eta T}{\mathcal{E}+{\mathcal{P}}_0}.$$

$$\frac{\zeta }{{\tau }_{\Pi }}$$

Full text

††journal: Computer Physics Communications

**PROGRAM SUMMARY

** Manuscript Title: (3+1)-dimensional dissipative relativistic fluid dynamics at non-zero net baryon density

Authors: Lipei Du, Ulrich Heinz

Program Title: BEShydro

Licensing provisions: GPLv3

Programming language: C++

Computer: Laptop, desktop, cluster

Operating system: GNU/Linux distributions, Mac OS X

Memory usage: For a $121\times 121\times 121$ grid, 0.57 GB (conserved variables), 0.003 GB (EOS4 tables), 0.035 GB (dynamical sources)

Keywords: Heavy-ion collisions, quark-gluon plasma, relativistic hydrodynamics, non-zero baryon density, baryon diffusion, Beam Energy Scan, BEShydro

Classification: 12 Gases and Fluids, 17 Nuclear Physics

External routines/libraries: GNU Scientific Library (GSL)

Nature of problem: (3+1)-dimensional dynamical evolution of hot and dense matter created in relativistic heavy-ion collisions using second-order dissipative relativistic fluid dynamics, including evolution of net baryon number and its dissipative diffusion current.

Solution method: Runge-Kutta Kurganov-Tadmor algorithm

Running time: A test run with baryon diffusion on a 2-dimensional grid $121\times 121$ uses 0.32 sec/time step, and on a 3-dimensional grid $121\times 121\times 121$ uses 46.08 sec/time step on a MacBook Pro with a 2.7 GHz Intel Core i5 processor and 8 GB 1867 MHz DDR3 memory.

1. 1 Introduction

2. 2 Relativistic hydrodynamics

3. 2.1 Equations of motion

4. 2.1.1 Conservation laws

5. 2.1.2 Evolution equations for the dissipative flows

6. 2.1.3 Evolution equations in Milne coordinates

7. 2.2 Transport coefficients

8. 2.2.1 Shear stress tensor

9. 2.2.2 Bulk viscous pressure

10. 2.2.3 Baryon diffusion current

11. 2.3 Equation of State (EoS)

12. 2.3.1 Construction of the EoS

13. 2.3.2 Equations of state implemented in BEShydro

14. 3 Numerical scheme

15. 3.1 The Kurganov-Tadmor algorithm

16. 3.2 Numerical derivatives

17. 3.3 Root finding with baryon current

18. 3.4 Regulation scheme

19. 4 Code validation with semi-analytical solutions

20. 4.1 The Riemann problem

21. 4.2 Bjorken flow

22. 4.3 Gubser flow

23. 4.3.1 Gubser coordinates

24. 4.3.2 Gubser flow with baryon diffusion

25. 4.4 Comparison to other codes

26. 4.5 Tests summary

27. 5 Baryon diffusion in an expanding QGP

28. 6 Conclusions and outlook

29. A Conservative form of the evolution equations

30. B Explicit form

1 Introduction

Over the last decade, second-order dissipative relativistic fluid dynamics (RFD) [1, 2, 3] has developed into a powerful and phenomenologically very successful tool for the description of the dynamical evolution of the hot and dense matter created in relativistic heavy-ion collisions [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. It is an effective theory for the macroscopic evolution of the conserved quantities of the system (energy, momentum, baryon number, strangeness, isospin and electric charge), coupled to dissipative flows whose dynamics is controlled by the competition between microscopic scattering (which moves the system closer to local thermodynamic equilibrium) and macroscopic expansion (which tends to drive the system away from local equilibrium). In heavy-ion collisions, dissipative fluid dynamics describes well the evolution of the quark-gluon plasma (QGP), a strongly coupled plasma that exhibits almost perfect liquid behavior [17, 18, 19, 7, 20, 21] but exists only at temperatures above about 150 MeV. To obtain quantitatively precise predictions for heavy-ion collisions, the fluid dynamic stage must be initialized with the output from some microscopic model for the pre-hydrodynamic stage at very early times, when the energy deposited in the collision fireball is still so far away from local momentum isotropy and thermal equilibrium that a hydrodynamic language can not yet be applied [22, 23, 24, 25, 26, 27], and it must be coupled to a microscopic kinetic description of the hadronic rescattering stage at the end of the fireball evolution [28, 29, 30, 31, 32], after the QGP has converted back into a gas of hadrons and hadronic resonances which quickly moves towards “kinetic freeze-out” when all strong interactions cease and the energies and momenta of all strongly interacting particles stop changing.

The initial development of dissipative RFD ignored the evolution of conserved currents such as net baryon number and strangeness because the community’s attention was focused on experiments performed at the highest available collision energies at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the Large Hadron Collider (LHC) at CERN, at which the colliding atomic nuclei are largely transparent to each other, creating a system of approximately zero net baryon number, strangeness, and isospin charge near midrapidity in the center of mass frame (formed by the interaction of low-$x$ gluons from the colliding nuclei which carry no conserved charges) whereas the net baryon number and electric charges brought into the collision by the valence quarks within the incoming nuclei cannot be stopped and end up mostly at far forward and backward rapidities where they are very hard to measure experimentally. Only with the Beam Energy Scan (BES) program at RHIC [33], in which heavy-ion collisions were studied at lower collision energies where some of the incoming baryon charge gets stopped near midrapidity, became the need urgent for including the dynamics of the baryon number and other conserved charge currents in the hydrodynamic description. These developments are also relevant for the theoretical description of future experiments at NICA [34] and FAIR [35, 36].

We here describe a new (3+1)-dimensional dissipative RFD code, which we call BEShydro, that includes the evolution of the net baryon charge and diffusion currents.111The code is open source and can be freely downloaded from https://github.com/LipeiDu/BEShydro BEShydro evolved out of the CPU version of the code GPU-VH [37], by adding a number of additional dissipative terms [16, 38, 15] describing the evolution of the bulk and shear viscous pressures, in addition to the evolution equations for the net baryon charge and diffusion currents. As far as we know, at this point in time the only other code that shares all of the main features of BEShydro is the latest version of MUSIC [39], while other codes (e.g. [14, 37, 40]) so far ignore the evolution of the net baryon diffusion current or, in some cases, even that of the net baryon charge. BEShydro has been developed completely independently of MUSIC; it can thus serve as a platform for detailed code validations and comparisons, even if in future applications the two codes will likely be applied to different collision systems, using different initialization modules and hadronic afterburners.

The purpose of this document is to review the physics ingredients needed for describing the hydrodynamic evolution of systems with non-zero conserved charges (such as baryon number) and to provide a detailed description of the structure and performance of BEShydro. We start from the most general form of the equations solved by BEShydro in Sec. 2, describe the numerical scheme used by the code in Sec. 3, and include a number of validation tests in Sec. 4, using specific, highly symmetric settings in which analytic or semi-analytic solutions of the hydrodynamic evolution equations have been found. Some of the novel tests of the charge transport sector of the code described here should be useful for the developers of other codes with the ability to describe dissipative effects related to charge diffusion. At this point BEShydro evolves only a single conserved charge (net baryon number); future generalization to include strangeness and isospin [41] is expected to be conceptually straightforward.

2 Relativistic hydrodynamics

In this section the physics ingredients of BEShydro are described: the propagated physical quantities, their equations of motion, and several variants of the Equation of State (EoS) controlling the expansion of the liquid formed in relativistic heavy-ion collisions.

2.1 Equations of motion

2.1.1 Conservation laws

Hydrodynamics is a macroscopic theory describing the space-time evolution of the 14 components of the energy-momentum tensor $T^{\mu\nu}(x)$ and the net (baryon) charge current $N^{\mu}(x)$.222For $N_{C}>1$ conserved charges, the number of evolved quantities increases to $10{\,+\,}4N_{C}$. Five evolution equations arise from the conservation laws for energy, momentum, and the baryon charge [15]:

[TABLE]

Here $d_{\mu}$ ($\mu=0,1,2,3$) stands for the covariant derivative in a general system of space-time coordinates, with metric tensor $g^{\mu\nu}$ defined with negative signature (“mostly minus” convention $(+,-,-,-)$), $g\equiv-\det{(g_{\mu\nu})}$, and the Christoffel symbols (see, e.g., [42])

[TABLE]

The 14 independent components of $T^{\mu\nu}$ and $N^{\mu}$ are more physically defined in terms of the hydrodynamic decomposition of these tensors [43],

[TABLE]

Here the flow 4-velocity $u^{\mu}(x)$, with $u^{\mu}u_{\mu}=1$, is defined as the time-like eigenvector of the energy-momentum tensor,

[TABLE]

and specifies the local rest frame (LRF) of the fluid at point $x$ (the so-called “Landau frame”). The tensors $u^{\mu}u^{\nu}$ and $\Delta^{\mu\nu}\equiv g^{\mu\nu}-u^{\mu}u^{\nu}$ then are projectors on the temporal and spatial directions in the LRF. ${\cal E}$ and ${\cal N}$ are the energy and net baryon density in the LRF which can be obtained as the following projections of $T^{\mu\nu}$ and $N^{\mu}$:

[TABLE]

From these, the local equilibrium pressure ${\cal P}_{0}$ is obtained through the EoS ${\cal P}_{0}={\cal P}_{0}({\cal E},{\cal N})$. The shear stress $\pi^{\mu\nu}$, the bulk viscous pressure $\Pi$, and the baryon diffusion current $n^{\mu}$ are dissipative flows describing deviations from local equilibrium.

Using the decomposition (4,5) the conservation laws (1,2) can be brought into the physically intuitive form [44]

[TABLE]

Here $D=u_{\mu}d^{\mu}$ denotes the time derivative in the LRF, $\theta=d_{\mu}u^{\mu}$ is the scalar expansion rate, $\nabla^{\mu}=\partial^{\langle\mu\rangle}$ (where generally $A^{\langle\mu\rangle}\equiv\Delta^{\mu\nu}A_{\nu}$) denotes the spatial gradient in the LRF, and $\sigma^{\mu\nu}=\nabla^{\langle\mu}u^{\nu\rangle}$ (where generally $B^{\langle\mu\nu\rangle}\equiv\Delta^{\mu\nu}_{\alpha\beta}B^{\alpha\beta}$, with the traceless spatial projector $\Delta^{\mu\nu}_{\alpha\beta}\equiv\frac{1}{2}(\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}+\Delta^{\nu}_{\alpha}\Delta^{\mu}_{\beta})-\frac{1}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta}$) is the shear flow tensor. While these equations clearly exhibit the physics in the LRF (which varies from point to point), BEShydro solves the conservation laws (1,2) in a fixed global computational frame. Their explicit form in the global frame is discussed in Sec. 2.1.3.

The 5 conservation laws (8-10) are sufficient to determine the energy and baryon density, ${\cal E}$ and ${\cal N}$, together with the 3 independent components of the flow velocity $u^{\mu}$, as long as the shear stress $\pi^{\mu\nu}$, the bulk viscous pressure $\Pi$, and the baryon diffusion current $n^{\mu}$ vanish.333Note that the shear stress is traceless, $\pi^{\mu}_{\mu}=0$, and both $\pi^{\mu\nu}$ and $n^{\mu}$ have only spatial components in the LRF, $u_{\mu}\pi^{\mu\nu}=\pi^{\mu\nu}u_{\nu}=u_{\mu}n^{\mu}=0$. $\pi^{\mu\nu}$, $\Pi$, and $n^{\mu}$ thus describe 5+1+3=9 dissipative degrees of freedom. Their evolution is not directly constrained by conservation laws but controlled by the competition between microscopic scattering processes (which drive the system towards local equilibrium and the dissipative flows to zero) and the macroscopic expansion (which drives the system away from equilibrium and the dissipative flows away from zero). Their evolution is thus controlled by both micro- and macroscopic physics. One way to obtain their evolution equations is DNMR theory [16, 38, 15] which uses the method of moments of the Boltzmann equation and which we employ here.

Expressed through its natural variables, i.e. the temperature $T$ and baryon chemical potential $\mu$, the equilibrium pressure ${\cal P}_{0}({\cal E},{\cal N})={\cal P}_{0}(T,\mu)$ is recognized as the grand-canonical thermodynamic potential for a system with temperature $T({\cal E},{\cal N})$ and chemical potential $\mu({\cal E},{\cal N})$. In principle, $T$ and $\mu$ are not needed for the hydrodynamic evolution, but they may be required to compute certain signatures of the evolving fluid (such as the spectrum of electromagnetic radiation emitted during its evolution or the spectrum of hadrons into which it decays at the end of the life of the quark-gluon plasma phase), and in BEShydro the driving force for net baryon number diffusion is formulated in terms of the gradient of $\mu/T$ rather than that of the net baryon density. Also, the transport coefficients controlling the evolution of the dissipative flows are most naturally expressed as functions of $T$ and $\mu$ since they are defined as response functions of the thermal equilibrium system described by the potential ${\cal P}_{0}(T,\mu)$. Different versions of the EoS ${\cal P}_{0}(T,\mu)$ or ${\cal P}_{0}({\cal E},{\cal N})$ used in BEShydro will be described in Sec. 2.3.

2.1.2 Evolution equations for the dissipative flows

In BEShydro the dissipative flows are evolved with DNMR theory [15, 38, 16]. While the equations of motion in this theory are derived from the Boltzmann equation which is applicable only to weakly coupled systems [45], the hydrodynamic description is an effective theory which is generic and applicable also in the strong coupling regime where the Boltzmann equation is not valid [46]. When applying the DNMR equations to the fluid produced in nuclear collisions, which appears to be strongly coupled, one must replace its material properties, i.e. the EoS and transport coefficients, by those for real QCD matter.

In the framework of DNMR theory, the dissipative transport equations are given by the following relaxation equations:

[TABLE]

Here $(Dn)^{\langle\mu\rangle}\equiv\Delta^{\mu\nu}Dn_{\nu}$ and $(D\pi)^{\langle\mu\nu\rangle}\equiv\Delta^{\mu\nu}_{\alpha\beta}D\pi^{\alpha\beta}$, ensuring that all terms are purely spatial in the LRF and, where applicable, traceless. $\tau_{\Pi}$, $\tau_{n}$, and $\tau_{\pi}$ are the relaxation times for $\Pi$, $n^{\mu}$, and $\pi^{\mu\nu}$, respectively. They control how fast the dissipative flows relax to their Navier-Stokes limits [44]:

[TABLE]

where $\zeta$, $\kappa_{n}$, and $\eta$ are the bulk viscosity, baryon diffusion coefficient, and shear viscosity, respectively, describing the first-order response of the dissipative flows to their driving forces, the (negative of the) scalar expansion rate $\theta$, the spatial gradient of $\mu/T$ in the LRF, $\nabla^{\mu}(\mu/T)$, and the shear flow tensor $\sigma^{\mu\nu}$, respectively, which drive the system away from local equilibrium.

The scalar, vector and tensor source terms on the r.h.s. of Eqs. (11-13), $\mathcal{J},\;\mathcal{K},\;\mathcal{R},\;\mathcal{J}^{\mu},\;\mathcal{K}^{\mu},\;\mathcal{R}^{\mu},\;\mathcal{J}^{\mu\nu},\;\mathcal{K}^{\mu\nu}$ and $\mathcal{R}^{\mu\nu}$ contain terms of second order in the small parameters Knudsen number (ratio between a characteristic microscopic and macroscopic time or length scale of the fluid) and inverse Reynolds number (ratio between dissipative quantities and local equilibrium values). According to the notation established in [16], the Navier-Stokes and $\mathcal{K}$ terms on the r.h.s. of Eqs. (11-13) are of first and second order in the Knudsen number(s), respectively, the $\mathcal{J}$ terms are of order Knudsen number times inverse Reynolds number, and the $\mathcal{R}$ terms are of second order in the inverse Reynolds number(s). Their explicit expressions can be found in Ref. [16]. Following the arguments in Ref. [37], we here include only a subset of the $\mathcal{J}$ terms. As BEShydro is ultimately designed for precision studies of relativistic heavy-ion collisions, one should perhaps not put too much blind trust into these arguments and rather check their validity; on the other hand, adding the missing second-order source terms to the code at a later time should be straightforward (even if additional code stability tests may be needed). Future code updates will include additional terms as required by specific applications.

As implemented in the code, the thus simplified relaxation equations for the dissipative flows read

[TABLE]

here $\omega^{\mu\nu}=\frac{1}{2}(\nabla^{\mu}u^{\nu}-\nabla^{\nu}u^{\mu})$ is the vorticity tensor. The additional transport coefficients $\delta_{\Pi\Pi}$, $\lambda_{\Pi\pi}$, $\tau_{n}$ etc. will be discussed in Sec. 2.2. In Eqs. (18) and (19) we removed the transversality constraints on the l.h.s. by using footnote 3 and

[TABLE]

as well as its analog for $\pi^{\mu\nu}$, moving the extra terms as additional source terms to the r.h.s.

It is worth pointing out that in Eqs. (17)-(19) we have followed Ref. [39] in ignoring terms describing the direct influence of baryon diffusion, $n^{\mu}$, on the evolution of the shear and bulk viscous stresses, $\pi^{\mu\nu}$ and $\Pi$. Baryon evolution still affects the evolution of the system indirectly through the EoS. In this approach it has been shown [39, 27] that, while dissipative baryon diffusion effects directly influence the net-proton distributions, its indirect effects on the distributions of mesons and charged hadrons are negligible. It might be interesting to study to which extent second-order couplings between baryon diffusion and viscous stresses can modify this conclusion.

2.1.3 Evolution equations in Milne coordinates

Up to this point the formalism is completely general. For application to ultra-relativistic heavy-ion collisions we need the specific form of the evolution equations in Milne coordinates $x^{\mu}=(\tau,x,y,\eta_{s})$ which are best adapted to the relativistic collision kinematics and subsequent almost boost-invariant longitudinal flow pattern [47]. In terms of Cartesian coordinates $(t,x,y,z)$ the longitudinal proper time $\tau$ and space-time rapidity $\eta_{s}$ are defined as

[TABLE]

The mid-rapidity point $z=\eta_{s}=0$ at $\tau=0$ defines the collision point in the global (computational) frame. In Milne coordinates the metric tensor is

[TABLE]

the fluid four-velocity is $u^{\mu}=(u^{\tau},u^{x},u^{y},u^{\eta})$, and the four-derivative is $\partial_{\mu}=(\partial_{\tau},\partial_{x},\partial_{y},\partial_{\eta})$.444In all sub- and superscripts $\eta$ is short for $\eta_{s}$. The metric has the following non-vanishing Christoffel symbols:

[TABLE]

Plugging them into Eqs. (1,2) we obtain the conservation laws in Milne coordinates:

[TABLE]

Introducing the convective time derivative $d\equiv u^{\mu}\partial_{\mu}$, the relaxation equations can be written as

[TABLE]

with the shorthand notations $G^{\mu\nu}_{\pi}\equiv u^{\alpha}\Gamma^{\mu}_{\alpha\beta}\pi^{\beta\nu}+u^{\alpha}\Gamma^{\nu}_{\alpha\beta}\pi^{\beta\mu}$ and $G^{\mu}_{n}=u^{\alpha}\Gamma^{\mu}_{\alpha\beta}n^{\beta}$ for the geometrical source terms obtained when splitting the covariant LRF time derivative $D$ in Eqs. (17)-(19) into the convective time derivatives $d$ and a remainder (for example, $Dn^{\mu}\equiv u^{\alpha}d_{\alpha}n^{\mu}=u^{\alpha}(\partial_{\alpha}n^{\mu}+\Gamma^{\mu}_{\alpha\beta}n^{\beta})$). The $I$-terms are explicitly

[TABLE]

with

[TABLE]

The conservation laws (24-28) together with the dissipative transport equations (29-31) constitute the equations of motion of the relativistic hydrodynamic system encoded in BEShydro. Next we will discuss the transport coefficients appearing in these equations.

2.2 Transport coefficients

The EoS and transport coefficients describe the medium properties of the expanding fluid and as such must be determined microscopically. While for the EoS detailed knowledge is available now from lattice QCD (see Sec. 2.3), the same is not true for the transport coefficients. We will here use rough estimates for the transport coefficients that have been obtained from kinetic theory, but have to leave their precise determination to future theoretical or phenomenological work. Specifically, BEShydro implements the transport coefficients from Ref. [16, 48] which starts from the Boltzmann equation in Relaxation Time Approximation (RTA) and employs the 14-moment approximation [38, 16] for a one-component gas of Boltzmann particles with non-zero but small mass $m\ll T$. In all expressions we only keep the leading non-zero term in powers of $m/T\ll 1$. This is motivated by the approximate masslessness of the microscopic quark-gluon degrees of freedom that make up the fluid described by BEShydro.

2.2.1 Shear stress tensor

For the evolution of the shear stress we take the transport coefficients

[TABLE]

Following [49] we express shear viscous effects in terms of the kinematic shear viscosity

[TABLE]

For $\mu=0$ the kinematic shear viscosity reduces to the specific shear viscosity $\eta/{\cal S}$, where ${\cal S}$ denotes the entropy density, but it differs greatly from $\eta/{\cal S}$ at large net baryon densities, and this can lead to significant differences in the hydrodynamic flow patterns [50]. Parametrizations of $\bar{\eta}$ as a function of $T$ and $\mu$ are discussed in [51, 52]; in BEShydro the default setting for $\bar{\eta}$ is a constant, $\bar{\eta}=0.2$. Given $\bar{\eta}$, Eq. (38) is used to calculate $\tau_{\pi}=5\bar{\eta}/T$, and Eqs. (39) to obtain the remaining transport coefficients.

2.2.2 Bulk viscous pressure

To evolve the bulk viscous pressure we take the transport coefficients

[TABLE]

where $c_{s}$ is the speed of sound in the medium (see Eq. (61)). Similar to the shear viscosity, we start from a parametrization as a function of $T$ and $\mu$ of the kinematic bulk viscosity

[TABLE]

and determine from it the bulk relaxation time $\tau_{\Pi}$ using Eq. (41) and the remaining transport coefficients using Eqs. (42). In the code, we use a parametrization that interpolates between lattice QCD data for the QGP phase and results obtained from the hadron resonance gas model for the hadronic phase, connected quadratically around the pseudocritical temperature $T_{c}=155$ MeV [53]:

[TABLE]

with $x=T/T_{c}$ and fitted parameters

[TABLE]

The bulk viscous pressure describes the deviation from the thermal pressure for a non-perfect expanding or contracting fluid. Bulk viscosity has been shown to generate important effects on the slope of the transverse momentum spectra and their azimuthal anisotropy [54]. Bulk viscous effects are expected to be strongest near the quark-hadron phase transition [55, 56, 57, 58, 59, 60, 61], especially near the QCD critical point where critically enhanced contributions associated with critical slowing-down play a key dynamical role [62, 63, 64].

2.2.3 Baryon diffusion current

Compared to the transport coefficients related to the bulk and shear stresses, those controlling baryon diffusion are much less explored. Following Ref. [39], we use the coefficients obtained from the Boltzmann equation for an almost massless ($m/T\ll 1$) classical gas (for a calculation of the baryon diffusion coefficient $\kappa_{n}$ for a massive gas of hadrons see Ref. [65]):

[TABLE]

Here $\tau_{n}$ is the relaxation time of the baryon diffusion current in Eq. (12) and parametrized as

[TABLE]

with a free parameter $C_{B}$.555Note that this procedure is different from the bulk and shear viscosity where we parametrized the first-order transport coefficients and computed the relaxation times from them, whereas here we parametrize the relaxation time and use it to compute the diffusion coefficient. The expression (45) for the diffusion coefficient $\kappa_{n}$ was derived in first-order Chapman-Enskog approximation [39] whereas the second-order transport coefficients (46) were obtained in the 14-moment approximation [48, 38, 16]. In the limit of small net baryon density ($\mu\to 0$) the diffusion coefficient $\kappa_{n}$ reduces to

[TABLE]

In RTA, where the collision term in the Boltzmann equation is parametrized with a single relaxation time $\tau_{\mathrm{rel}}$, the relaxation times for $\Pi$, $\pi^{\mu\nu}$ and $n^{\mu}$ are all the same, i.e. $\tau_{\Pi}=\tau_{\pi}=\tau_{n}=\tau_{\mathrm{rel}}$. Here we allow them to be different. If one does, however, impose the constraint $\tau_{n}=\tau_{\pi}=\tau_{\mathrm{rel}}$, the baryon diffusion and shear viscosity coefficients can be related as follows [66]:

[TABLE]

In contrast to Eqs. (45,48), this expression takes into account quantum statistics. The function $C(\mu/T)$ exhibits a weak dependence on $\mu/T$, interpolating between 5/3 at large $\mu/T$ and a somewhat smaller value at small $\mu/T$ whose precise magnitude depends on the number of massless degrees of freedom in the gas [66]. Note that, since for small $\mu/T$ the net baryon density ${\cal N}\propto\mu$, both (48) and (49) yield nonzero baryon diffusion coefficients at zero net baryon density. At large $\mu/T$, the authors of Ref. [66] have shown that the ratio $\kappa_{n}T/\eta$ approaches zero, i.e., at large net baryon densities and low temperatures baryon diffusion effects can generally be neglected in comparison with shear viscous stresses.

In addition to kinetic theory, gauge/gravity duality has also been widely used to determine the transport properties of the QGP (see, for example, [67, 68, 69]). With this method it is possible to study the transport properties of strongly coupled gauge theories for which no kinetic theory description exists, including the critical dynamics near a critical point. The latter is expected to provide critical signatures for the experimental identification of the QCD critical point [62, 64]. Studying and comparing the baryon number evolution with baryon diffusion coefficients corresponding to a weakly [39] or a strongly coupled QGP [68] can be interesting [27, 70].

2.3 Equation of State (EoS)

Another important medium property that crucially affects the dynamical evolution of the fluid is its equation of state

[TABLE]

In practice, for the calculation of the transport coefficients and chemical forces we also need the equivalent relations $T({\cal E},{\cal N})$ and $\mu({\cal E},{\cal N})$. We use the term EoS generically for any one of these relations.

2.3.1 Construction of the EoS

Since the matter produced in nuclear collisions passes through very different physical regimes that differ by orders of magnitude in energy density and must be described with different effective degrees of freedom, we need an EoS that describes the medium properties over a wide range of temperature and length scales [71]. On most scales the degrees of freedom of the evolving system are strongly coupled, rendering perturbative investigations from first principles unreliable. Over the last decades, lattice QCD (LQCD) has been established as the most precise non-perturbative framework to calculate the EoS of strongly interacting matter at zero baryon chemical potential (see, e.g., [72]). The method works well at temperatures above $\sim 100$ MeV; at lower temperatures the lattice signals become weaker and more noisy, necessitating the matching of LQCD data to an analytical hadron resonance gas model.

Lattice QCD obtains the EoS by calculating the trace of the energy-momentum tensor $T^{\mu\nu}$, ${\cal E}-3{\cal P}_{0}$ (usually referred to as the “trace anomaly” or “interaction measure”), describing deviations from the conformal EoS. Defining the rescaled dimensionless trace anomaly

[TABLE]

the thermal pressure at zero chemical potential can be written as

[TABLE]

Unfortunately, this method cannot be directly extended to non-zero chemical potential where the evaluation of the QCD path integral for the interaction measure $I(T,\mu)$ suffers from a “sign problem” [71], precluding its direct computation with standard Monte-Carlo methods. This problem can be partially circumvented by using standard LQCD techniques to also compute the $\mu$-derivatives of the pressure ${\cal P}_{0}(T,\mu)$ at $\mu=0$ and construct ${\cal P}_{0}(T,\mu)$ at non-zero $\mu$ from its Taylor series around $\mu=0$:

[TABLE]

The expansion coefficients

[TABLE]

are known as the “baryon number susceptibilities” [73, 74, 75]. The computational effort of computing them increases rapidly with their order $n$; at this time, the Taylor expansion (53) includes terms up to order $n_{\mathrm{max}}=3$ and, near $T=T_{c}$, converges well up to about baryon chemical potentials $\mu/T\lesssim 2$ [73, 74, 75].

At low temperatures $T\ll T_{c}\simeq 155$ MeV, LQCD is increasingly affected by lattice artifacts and the system is more properly described in terms of hadronic degrees of freedom as a “hadron resonance gas” (HRG). In the HRG model, the interactions among different hadronic species are accounted for by including all experimentally identified scattering resonances as additional, non-interacting particle species. In the HRG model the interaction measure $({\cal E}{-}3{\cal P}_{0})/T^{4}$ is given as [76, 77]

[TABLE]

where particle species with spin-isospin degeneracy $g_{i}$ and mass $m_{i}$ smaller than some cut-off $m_{\mathrm{max}}$ can be included, and $\eta_{i}=-1~{}(+1)$ for bosons (fermions) describes the effects of quantum statistics.

In principle, the EoS used for the hydrodynamic evolution should include the same set of hadronic resonances as the hadronic afterburner employed to describe the kinetic final freeze-out stage because otherwise a mismatch of the energy and baryon densities occurs on the conversion surface where we change between these two different dynamical descriptions. In practice these discontinuities tend to be small, and such care is not always taken. In our applications of BEShydro we use different matched equations of state for different hadronic afterburners (e.g., for UrQMD [28, 29] and SMASH [30]); however, the module for matching the lattice QCD data to a HRG with adjustable hadronic mass spectrum is not part of this code distribution.

For the matching procedure between the LQCD and HRG equations of state different methods have been used. For example, in Refs. [74, 39] the pressure is interpolated as follows:

[TABLE]

Here ${\cal P}_{0}^{\mathrm{HRG}}$ and ${\cal P}_{0}^{\mathrm{LAT}}$ are the equilibrium pressures for the hadron resonance gas and from lattice QCD, respectively, $T^{\prime}(\mu)$ is the “switching temperature” and $\Delta T^{\prime}$ controls the width of the “overlap region”. The authors of Ref. [78], on the other hand, interpolate the interaction measure at $\mu=0$, $I(T,\mu=0)$ smoothly between $T_{1}=T_{c}(\mu{=}0)=155$ MeV and $T_{2}=180$ MeV, using a polynomial interpolation function.

Once the pressure ${\cal P}_{0}(T,\mu)$ is given, other thermodynamic quantities can be calculated from thermodynamic identities:

[TABLE]

The last equation requires expressing ${\cal E}$ and ${\cal N}$ through $T$ and $\mu$ after taking the derivatives. In practice, the functions ${\cal E}(T,\mu)$ and ${\cal N}(T,\mu)$ are numerically inverted, and the quantities $T,\mu,{\cal S},{\cal P}_{0}$, and $c^{2}_{s}$, as well as the two derivatives on the r.h.s. of Eq. (61), are stored in a table on a $({\cal E},{\cal N})$ grid which is interpolated by the hydrodynamic code as needed.

2.3.2 Equations of state implemented in BEShydro

In BEShydro four different equations of state, EOS1 to EOS4, are implemented, for different purposes: at zero chemical potential, we include a massless (conformal) EoS (${\cal E}=3{\cal P}_{0}$, EOS1) as well as an interpolated LQCD-HRG EoS from the Wuppertal-Budapest collaboration [72] (EOS2); at non-zero chemical potential, an appropriately generalized conformal EoS (EOS3) and an interpolated LQCD-HRG EoS from Ref. [39] (EOS4) are used.

EOS1 assumes an ideal gas of massless quarks and gluons:

[TABLE]

where $N_{c}=3$ and $N_{f}=2.5$ are the numbers of colors and (approximately) massless quark flavors, respectively.666Strange quarks, whose mass is of the same order of magnitude as the quark-hadron transition temperature, are (somewhat roughly) counted as 1/2 massless quark flavor. In EOS1, $c_{s}^{2}=1/3$ and $\mu=0$. While the conformal EoS does not properly describe the properties of the matter produced in nuclear collisions, it is, owing to its simplicity, very useful for code testing. Technically, EOS1 can be used in BEShydro even when the baryon density and baryon diffusion currents are being evolved; in that case, these currents do not couple to the rest of the hydrodynamic system and evolve purely as background fields.

EOS2 from the Wuppertal-Budapest collaboration [72] can be used for realistic simulations at vanishing $\mu$, i.e. for heavy-ion collisions at ultra-relativistic collision energies with $\sqrt{s_{\mathrm{NN}}}\gg 100$ GeV, especially near mid-rapidity. More details about EOS2 can be found in Refs. [72, 37].

EOS3 is the generalization of EOS1 to non-zero $\mu$. Starting from the ideal massless parton gas expression [79]

[TABLE]

where $N_{c}=3$ and the sum goes over massless quark flavors, we simplify it by setting $\mu_{f}=\mu/3$ for all flavors (which is appropriate if only baryon number is considered as a conserved charge):

[TABLE]

with $p_{0}=(16+10.5N_{f})\pi^{2}/90$. We again count strange quarks with a factor 1/2, i.e. we set $N_{f}=2.5$ so that for $\mu=0$ (64) reduces to (62). Inverting these functions numerically one obtains the EoS tables used in the hydrodynamic code (see Fig. 1).

EOS4 from Ref. [39], extended to finite baryon chemical potential by combining a lattice EoS at high temperature and a HRG EoS at low temperature with a Taylor expansion in $\mu/T$ using techniques discussed in Sec. 2.3.1, allows to study the evolution of systems with non-zero net baryon density. It is plotted in Fig. 2. In the code, tabulated values for ${\cal P}_{0},\,T$ and $\mu/T$ as functions of $({\cal E},{\cal N})$ are included. If the code requires the EoS at $({\cal E},{\cal N})$, ${\cal P}_{0},T$ and $\mu/T$ are calculated from nearest neighbors in the table using 2D bilinear interpolation. We note that EOS4 does not include a critical point or first-order phase transition at large $\mu$. A lattice QCD based EoS that includes these features, with adjustable location of the critical point and strength of the first-order transition beyond that point, was constructed by the BEST Collaboration [74] and could be imported into BEShydro for future dynamical simulations aiming at helping to locate the QCD critical point.

3 Numerical scheme

We now describe the numerical scheme used in BEShydro to solve the coupled set of evolution equations (the conservation laws (24-28) together with the dissipative relaxation equations (29-31) and the EoS (50)) discussed in Section 2. Initial values for all components of the baryon charge current and energy-momentum tensor are set on a surface of constant longitudinal proper time $\tau$.777A dynamical initialization routine with sources for the divergences of the baryon current and energy-momentum tensor that describe the gradual “hydrodynamization” of the matter produced in the collision [26, 80, 27] will be discussed elsewhere. We focus our attention on aspects of the algorithm related to the evolution and influence of the baryon density and diffusion currents, referring interested readers to Refs. [8, 15, 37] for additional technical details.

3.1 The Kurganov-Tadmor algorithm

Using the definition of $d\equiv u^{\mu}\partial_{\mu}$, the equations of motion in Eqs. (24-28), (29-31) can all be written [15, 37] in the same first-order flux-conserving form

[TABLE]

where $v^{i}\equiv u^{i}/u^{\tau}$ ($i=x,y,\eta_{s}$) is the 3-velocity of the fluid, the conserved quantity $q$ can be any component (or linear combination of components) of $T^{\mu\nu}$ and $N^{\mu}$, and $S_{q}$ is the corresponding source term. The explicit equations (66) implemented in BEShydro are given in Appendix A. This form allows all quantities to be evolved with the same numerical transport algorithm.

BEShydro is designed for flexibility so that different physical limits can be easily studied. Only the propagation of $T^{\tau\tau}$, $T^{\tau x}$, $T^{\tau y}$, and $T^{\tau\eta}$ is always turned on. The evolution of all other variables (dissipation and/or baryon related) can be conveniently turned on and off independently; only the propagation of $n^{\mu}$ requires the evolution of $N^{\tau}$ to be turned on.

To solve equation (66) we use the Kurganov-Tadmor (KT) algorithm [81], with a second-order explicit Runge-Kutta (RK) ordinary differential equation solver [82] for the time integration step. This scheme is widely used in relativistic hydrodynamic simulations (see, e.g., [8, 40, 37]).

3.2 Numerical derivatives

For the source terms $S_{q}$ we must evaluate spatial and temporal derivatives of the thermodynamic variables and dissipative flows. For the time derivatives the code uses first-order forward differences:

[TABLE]

Here $A$ is the quantity to be differentiated, $i,j,k$ are integer labels for the $x,y$, and $\eta_{s}$ coordinates of the grid point, and $\Delta\tau$ is the temporal grid size (numerical resolution in the $\tau$ coordinate). $n$ and $n-1$ are temporal indices denoting the present and preceding time step. To initialize the temporal evolution code at the first time step $n=1$ we set $A^{0}_{i,j,k}=A^{1}_{i,j,k}$. This is especially important when the initial flow velocity $u^{\mu}$ is non-zero, for example in the case of the Gubser flow test described in Sec. 4.3.

The code provides two methods for calculating spatial derivatives. The first uses second-order central differences, i.e. the derivative of any quantity $A$, say, with respect to $x$ is calculated as

[TABLE]

where $\Delta x$ is the numerical resolution (grid size) in $x$ direction. The boundary conditions are taken care of by introducing ghost cells on the boundary as described in [37].

The second method calculates the spatial derivative from a combination of second-order central and first-order backward and forward derivatives, using the generalized minmod flux limiter:

[TABLE]

where the multivariate minmod function is defined as

[TABLE]

with $\mathrm{\texttt{minmod}}(x,y)\equiv[\mathrm{\texttt{sgn}}(x)+\mathrm{\texttt{sgn}}(y)]\cdot\mathrm{\texttt{min}}(|x|,|y|)/2$ and $\mathrm{\texttt{sgn}}(x)\equiv|x|/x$. In other words, $\mathrm{\texttt{minmod}}(x,y,z)$ always gives the value which is the closest to 0 among $(x,y,z)$. The parameter $\theta_{\mathrm{f}}\in[1,2]$; $\theta_{\mathrm{f}}=1$ ($\theta_{\mathrm{f}}=2$) corresponds to the most (least) dissipative limiter. In BEShydro Eq. (69) is used only for the derivatives of $u^{\mu}$ and ${\cal P}_{0}$, and only when selected by the user as an option.

3.3 Root finding with baryon current

The code evolves the components of the energy momentum tensor in the global computational frame, but the EoS (which is needed to close the set of evolution equations) and the computation of the source terms on the r.h.s. of Eq. (66) require knowledge of fluid velocity $u^{\mu}$ and the energy and baryon density in the local rest frame of the fluid. Computing the latter from the former is known as the “root finding” problem. This must be done as efficiently as possible since this problem must be solved at every point of the computational space-time grid.

At finite baryon density, with nonzero baryon diffusion current, the root finding algorithm becomes more complex than described in Ref. [37]. We here describe the most general form of the root finding problem: assuming that $T^{\tau\mu}$, $N^{\tau}$, $\pi^{\tau\mu}$, $\Pi$, and $n^{\tau}$ are all known from the latest temporal update step, we want to compute ${\cal E}$, ${\cal N}$, and $u^{\mu}$. As will be demonstrated in Sec. 4, the following algorithm [14, 83, 40] works for both ideal and dissipative fluids, i.e. for both vanishing and non-vanishing dissipative flows. We start by introducing the “ideal fluid contributions” $M^{\mu}$ and $J^{\tau}$ to the energy-momentum current $T^{\tau\mu}$ and baryon density $N^{\tau}$ in the computational frame:

[TABLE]

Note that for a viscous fluid ${\cal P}={\cal P}_{0}+\Pi$ includes implicitly the bulk viscous pressure. We use the following only when baryon evolution is turned on; otherwise, we use the simpler algorithm described in [37] where ${\cal E}$ is found first, using a 1-dimensional zero search. Here we first find the magnitude of the flow velocity, $v$, by solving iteratively [14, 83, 40]

[TABLE]

where $M\equiv\sqrt{(M^{x})^{2}+(M^{y})^{2}+\tau^{2}(M^{\eta})^{2}}$ and ${\cal E}(v)$, ${\cal N}(v)$ are obtained from the known quantities $M^{\tau}$, $M$, and $J^{\tau}$ as

[TABLE]

Once the flow magnitude $v$ is known, we also know the flow 4-velocity:

[TABLE]

for $i=x,y,\eta_{s}$. Note that the algorithm makes active use of all the numerically known components listed above.

Rewriting Eq. (74) in the form

[TABLE]

we solve it by the standard Newton-Raphson method, by repeatedly updating the velocity with

[TABLE]

where

[TABLE]

until a sufficiently accurate value is reached. The last term in (81) is evaluated numerically using

[TABLE]

where (for equations of state like EOS4) the derivatives $\partial{\cal P}_{0}/\partial{\cal E}$ and $\partial{\cal P}_{0}/\partial{\cal N}$ must be interpolated from the values stored in the EoS table to the pair $({\cal E},{\cal N})$ tried in each step of the iteration.

When $v$ gets close to the speed of light, the Newton-Raphson iteration must be modified to avoid excursions into the causally forbidden region $v>1$. This can lead to numerical instabilities and/or poor convergence. We therefore follow the recipe proposed in [83] and use Eq. (79) to solve for $v$ only if in the previous time step, at the spatial grid point in question, $v\leq 0.563624$ or, equivalently (see (77)), $u^{\tau}\leq 1.21061$. Otherwise we instead solve for $u^{\tau}$ (which has no upper limit), by employing the Newton-Raphson algorithm to find the zero of

[TABLE]

In this case we need to evaluate in each iteration

[TABLE]

where

[TABLE]

with $v=\sqrt{1-1/(u^{\tau})^{2}}$. In Sec. 4.3 it will be shown that the switch between the two schemes, as implemented in BEShydro, works seamlessly, moving smoothly from $v=0.563624$ to $u^{\tau}=1.21061$ across the switching point.

When interpolating the EoS table to obtain the derivatives needed on the r.h.s. of Eqs. (82,85) one can encounter numerical errors in regions where derivatives of the EoS change discontinuously, e.g. in the recently developed BEST EoS [74] which adds (using a certain prescription) a critical point and first-order phase transition to the LQCD-HRG interpolated EOS4. For such situations, BEShydro offers another option for the root finding that avoids calculating these derivatives, at the price of somewhat degraded convergence which can slow the root finding algorithm compared to the Newton-Raphson method. This second regulation scheme may also be preferred when evolving more than one conserved charge, in which case not having to compute the thermodynamic derivatives may overcompensate for the slower intrinsic convergence of the root finding algorithm.

The modified root finder employs the following simple iteration scheme: starting with an initial guess $v_{i}$ for the velocity (e.g. the solution at this grid point from the preceding time step), we determine \bigl{(}{\cal E}(v_{i}),{\cal N}(v_{i})\bigr{)} from Eqs. (75,76) and the EoS ${\cal P}_{0}({\cal E},{\cal N})$, compute an updated value $v_{i+1}$ of the velocity from

[TABLE]

and iterate these steps until convergence is reached. For $v\geq 0.563624$ or $u^{\tau}\geq 1.21061$, one instead updates $u^{\tau}$ using the equation

[TABLE]

until convergence is reached.

In principle, the two methods are equivalent and should find the same root, within the prescribed numerical precision. In the Gubser test described in Sec. 4.3 they are indeed shown to yield identical numerical results. In both methods, higher numerical precision of the solution should be demanded when solving for $v$, due to the speed limit $v<1$.

We point out that extra hydrodynamic variables are propagated in the code that are not used in the root finding algorithm, such as $\pi^{xx}$, $\pi^{xy}$, $\pi^{x\eta}$, and $n^{\eta}$. In principle, these could be computed from the other components of the shear stress and baryon diffusion current by using the tracelessness of $\pi^{\mu\nu}$ and the orthogonality of $\pi^{\mu\nu}$ and $n^{\mu}$ to $u^{\mu}$. Instead, we propagate all shear stress and baryon diffusion components dynamically and use the tracelessness and orthogonality conditions to check the numerical precision of the code.

3.4 Regulation scheme

The solution of the hydrodynamic equations of motion on discretized grids, the bilinear interpolation of the EoS table at each grid point, the iterative nature of the root finding algorithm, and the need for calculating derivatives numerically all engender unavoidable numerical errors. In addition, the numerical solution for the hydrodynamic variables can make excursions into regions where the approximations under which the evolution equations were derived (such as ignoring higher-order gradient terms) are no longer valid, and the numerical evolution algorithm produces unphysical results. This happens, in particular, because Nature provides us with initial conditions that exhibit unavoidable quantum fluctuations which can lead to local excursions outside the region of validity of dissipative hydrodynamics. Although dissipation usually erases such large fluctuations over short time scales [83, 37], this may not happen quickly enough to avoid numerical instability of the evolution algorithm.888Note that we are not even talking about thermal fluctuations during the hydrodynamic evolution (our code solves deterministic equations of motion) which add possibly large stochastic fluctuations throughout the evolution history [84, 85].

In realistic event-by-event simulations such fluctuations can result in large gradients in both space and time which the code has to be able to deal with. Large gradients of the macroscopic variables can yield locally large Knudsen numbers (for which the fluid dynamic approximation breaks down) or large inverse Reynolds numbers (in which case the applicability of the 14-moment approximation used to simplify the hydrodynamic equations of motion is doubtful) [83, 37]. To ensure numerical stability of the code, such excursions must be regulated. To avoid the undesirable consequence that, after regulation, the algorithm no longer solves the underlying evolution equations, the regulation must be local, i.e. it must affect only very localized space-time regions, and its effects must be monitored so that the user is warned when the regulation becomes so strong and the regulated regions become so large that the code no longer correctly simulates the physics encoded in the evolution equations.

In practice, large gradients can drive large shear stress, bulk viscous pressure and baryon diffusion currents, and these can result in numerical instability or failure of the root finding algorithm. When this happens it is typically during the earliest evolution stage (where both the physical inhomogeneities driven by quantum fluctuation and the longitudinal expansion rate are largest) and/or in the very dilute regions near the transverse and longitudinal edge of the computing grid where the dissipative corrections to the leading thermodynamic quantities become large and the matter can no longer be reasonably treated as a fluid. Since the latter regions are typically far outside the domain where the matter is in the quark-gluon plasma phase (and thus outside the region where we want to apply the hydrodynamic picture), regulating them is innocuous as long as the regulation effects do not have sufficient time to propagate back inwards into the QGP region. Regulating large initial fluctuations is more dangerous because the fluctuations can be large both in- and outside the QGP phase. Both types of regulations must be carefully monitored.

Regulation schemes can be tricky, and a variety of implementations exist.999For example, CLVisc [40] requires $\texttt{max}(|\pi^{\mu\nu}|)<T_{0}^{\tau\tau}$; when this is violated for some cell in the dilute region, $\pi^{\mu\nu}$ is set to 0 locally. vHLLE [14] requires $\texttt{max}(|\pi^{\mu\nu}|/|T_{0}^{\mu\nu}|)<C$ and $|\Pi|/{\cal P}_{0}<C$, with $C$ being a constant of order but smaller than 1; if one of these conditions is violated, $\Pi$ and/or $\pi^{\mu\nu}$ are rescaled by a factor (which is common for all components of $\pi^{\mu\nu}$) to satisfy this requirement. In BEShydro, we follow the lead of iEBE-VISHNU [83], GPU-VH [37] and MUSIC [39] and implement two types of regulation that build on the schemes suggested in these earlier codes. Both are triggered by large dissipative flows which are then regulated. The trigger criterion compares (in ways defined more precisely below) $\pi^{\mu\nu}$ with $T_{0}^{\mu\nu}={\cal E}u^{\mu}u^{\nu}-{\cal P}_{0}\Delta^{\mu\nu}$, $\Pi$ with $\sqrt{{\cal E}^{2}+3{\cal P}_{0}^{2}}$, and $n^{\mu}$ with ${\cal N}u^{\mu}$.

For the shear stress tensor, iEBE-VISHNU [83] and GPU-VH [37] require

[TABLE]

with $\rho_{\mathrm{max}}\leq 1$. In addition, the tracelessness of $\pi^{\mu\nu}$ and its orthogonality to $u_{\nu}$ are required,

[TABLE]

where $\sqrt{\pi^{\mu\nu}\pi_{\mu\nu}}$ sets the scale and $\xi_{0}\ll 1$ is a small number [83]. At grid points where these trigger conditions are violated, Refs. [83, 37] regulate the shear stress tensor $\pi^{\mu\nu}$ by (see the left plot of Fig. 3)

[TABLE]

where

[TABLE]

For the bulk viscous pressure, which can make the root finding process fail when negative and too large, Refs. [83, 37] regulate $\Pi$ during the root finding process to ensure existence of at least one non-negative solution for ${\cal E}$, $v$ or $u^{\tau}$. In BEShydro we demand instead that

[TABLE]

and where this trigger condition is violated we regulate $\Pi$ by

[TABLE]

with

[TABLE]

In the same spirit, we require for the (space-like) baryon diffusion current

[TABLE]

In the code we replace these conditions by

[TABLE]

When one of these conditions is violated in a cell it triggers the following regulation of the baryon diffusion current:

[TABLE]

with

[TABLE]

Equations (90,93,97) define the first of our two regulation schemes. In contrast to Refs. [83, 37] where $\Pi$ was regulated during the root finding process, our regulation here is performed only after finishing each step of the two-step RK-KT algorithm. After implementing it we find that there is no need for additional regulation of $\Pi$ during the root finding process. As defaults for the regulation parameters we use the same values $\xi_{0}=0.1$ and $\rho_{\mathrm{max}}=1$ as proposed in [83, 37]. For documentation of the sensitivity studies leading to these default values we refer the reader to Ref. [83].

A second, different regulation scheme is based on the one implemented in the latest version of MUSIC [39]. (A flag in the code allows the user to select the preferred regulation scheme before running it.) In this second scheme, the dissipative components are rescaled by

[TABLE]

where regulation of quantity $i$ is triggered whenever $r_{i}$ exceeds the corresponding maximally allowed value $r_{i}^{\mathrm{max}}$ $(i=\pi,\Pi,n)$,101010If the regulation for the shear stress or baryon diffusion is triggered, all components of $\pi^{\mu\nu}$ or $n^{\mu}$ are regulated by a common regulation factor. with $r_{i}$ defined by

[TABLE]

By comparing Eq. (99) to Eqs. (90,93,97) and Eq. (100) to Eqs. (91,94,98) above, we see that the factors $r_{i}^{\mathrm{max}}/r_{i}$ play the same role as $\tanh{\rho_{i}}/\rho_{i}$ in the first regulation scheme, causing stronger regulation for larger $r_{i}$, with $r_{i}$ playing the role of the quantity $\rho_{i}$ whereas $f_{s}$, defined by

[TABLE]

playing a similar role as $\rho_{\mathrm{max}}$: as $f_{s}$ or $\rho_{\mathrm{max}}$ grows larger, the regulation gets weaker. $f_{s}$ is designed to approach $\chi_{0}$ when ${\cal E}\gg{\cal E}_{0}$ and 0 when ${\cal E}\ll{\cal E}_{0}$.

The right plot of Fig. 3 shows that for ${\cal E}<{\cal E}_{0}$, i.e. in the dilute region, $f_{s}$ decreases exponentially and the regulation strength increases accordingly. On the other hand, for large values of the parameter $\chi_{0}$, the regulation will hardly ever be triggered in the dense region ${\cal E}>{\cal E}_{0}$. Thus, unlike the first method, which always regulates larger dissipative components more strongly, irrespective of the energy density at the grid point, for the choice ${\cal E}_{0}=0.1$ GeV/fm3 [39] the second method causes hardly any regulation at grid points in the dense QGP region but more frequent and stronger regulation in the dilute region far outside the QGP fluid. The authors of Ref. [39] used $\chi_{0}=10$, and chose $r^{\mathrm{max}}_{n}=1$ for regulating the baryon diffusion current; we adopt the same value for $\chi_{0}$ and identical maximum $r$ values for all dissipative flows: $r^{\mathrm{max}}_{\pi}=r^{\mathrm{max}}_{\Pi}=r^{\mathrm{max}}_{n}=1$.

The default set for the regulation parameters is not universal and may need adjustments for different initial conditions, collision systems, and collision energies. The user is encouraged to play with these parameters to achieve maximal code stability with minimal changes to the physics encoded in the evolution equations. The regulation scheme may need to become more involved in future versions of dissipative hydrodynamics that include possibly large thermal and/or critical fluctuations in the dynamics (see e.g. [84, 85]).

4 Code validation with semi-analytical solutions

Numerical codes solving second-order (“causal”) relativistic dissipative fluid dynamics in 3+1 dimensions have only been developed over the last decade. They solve a problem for which in general no analytic solutions are available. Careful validation of any such code by testing its various components in simplified settings for which analytic or semi-analytic solutions are available is therefore mandatory. Some of these tests are nowadays standard and are included with this distribution precisely for that reason. The precursor of this code, the CPU-version of GPU-VH [37] was carefully validated using similar tests, but the entire baryon evolution sector in BEShydro is new such that direct comparisons with GPU-VH are of limited value. We therefore include here in particular novel semi-analytic tests of the baryon evolution equations. Direct code-to-code comparisons with MUSIC (whose latest version [39] also includes baryon evolution) will become possible when that version of MUSIC becomes public.

Building on validation protocols described in Refs. [14, 83, 37, 40, 39], we here discuss tests in which we compare, for identical initial conditions, the numerical solutions from BEShydro with (semi-)analytic solutions using Mathematica [86] for the Riemann problem for the Euler equations [87, 88, 89, 90, 91], for Bjorken flow [47], and for Gubser flow [92, 93] extended to systems with non-zero net baryon density and baryon diffusion current induced by a fluctuation in the initial state. We also include a direct comparison of BEShydro with the independently developed numerical algorithm described in Ref. [94] for a system with non-zero net baryon density in a (1+1)-dimensional setting with general longitudinal but vanishing transverse flow. By generalizing previously developed validation protocols to systems with non-zero net baryon density and baryon diffusion currents the work described in this section prepares the ground for code validation of other hydrodynamic codes at finite baryon density that are presently being developed elsewhere for the study of heavy-ion collisions at BES energies.

All tests described in this section are done without code regulation, i.e. all the regulation schemes described in Sec. 3.4 are turned off. In the code, all dimensionful quantities are represented by numbers given in length units, using the appropriate powers of [fm]; when plotting the results we sometimes convert them to physical units by multiplying with the appropriate powers of $\hbar c=0.197$ GeV fm.

In the tests and all other applications of the code completed to date we have used the following grid spacings: $\Delta x=\Delta y=0.05$ fm and $\Delta\eta_{s}=0.02$; $\Delta\tau$ is adjusted as needed and described in each case below. In realistic simulations the choice of grid spacing has to be a compromise between computational economy and capturing relevant physical information (e.g. large gradients, especially at early times). For smooth, ensemble-averaged initial conditions, larger grid spacings than those listed above can be used, whereas simulating small collision systems (such as proton-proton collisions) may require finer grids.

4.1 The Riemann problem

We start with testing the code against an analytical solution of the Riemann problem for the Euler equations, which historically has played an important role in fluid dynamics research and in the development of hydrodynamic codes [87, 89]. Here a special case of the relativistic Riemann problem, known as Sod’s shock tube problem [88, 89, 90, 91], is considered describing the 1D evolution of an ideal fluid in the transverse plane during the decay of a discontinuity across a $(y,z)$ surface placed at $x{\,=\,}0$, separating two constant initial states (“left” ($l$) for $x<0$ and “right” ($r$) for $x>0$) at rest:

[TABLE]

In the longitudinal ($z$) direction the fluid is assumed to expand with a boost-invariant velocity profile $u^{z}/u^{\tau}{\,=\,}z/t$ (i.e. $u^{\eta}{\,=\,}0$), and the system is initialized along a surface of constant longitudinal proper time $\tau_{0}$.

The decay of this initial discontinuity gives rise to general features of the Riemann problem, characterized by three elementary waves. Two of them are a rarefaction wave and a shock, moving into the two initial state regions of high and low density, respectively. Between them, two additional states emerge, separated by the third wave, which is a contact discontinuity moving with the fluid [87, 89] (see Fig. 4(b)). For a conformal EoS, an analytical solution for this problem can be derived from the conservation laws with the boundary condition that across the contact discontinuity pressure and velocity are constant (Figs. 4(a,c)) while the density has a jump (Fig. 4(b)) [89, 15, 95]. For non-conformal equations of state at non-zero net baryon density no general analytic solution is known, and the numerical solution can give rise to complex features (see, e.g., [89, 96]). Since for EOS1 $\partial{\cal P}_{0}/\partial{\cal N}$ (which is needed in Eqs. (82,85) for the velocity finding algorithm) is zero, the evolution of baryon density decouples from that of the energy density and pressure.

In the code, the assumption of longitudinal boost-invariance is implemented by setting the number of cells in the longitudinal direction to 1. The initial profiles of ${\cal E}$ and $T$ are obtained from the EoS. For ideal fluid dynamics the equations of motion become

[TABLE]

They are solved with initial conditions (102) with the default parameters listed in Fig. 4. For details about the analytical solution we refer the reader to Ref. [95] (see also [90, 91]). We point out that in the analytical solution from Ref. [95] the evolution starts at time zero whereas in the code the hydrodynamic evolution is initialized at $\tau_{0}=0.5$ fm. Since the solution is self-similar and depends only on the variable $x/(\tau-\tau_{0})$, we therefore compare in Fig. 4 the numerical results at $\tau=8.5$ fm/$c$ to the analytical solution at $\tau-\tau_{0}=8.0$ fm/$c$. For simplicity, the numerical test is done in Cartesian coordinates where all Christoffel symbols vanish. Also, to adequately capture the large discontinuity in the initial state, derivatives should be evaluated using Eq. (68), and not Eq. (69) which would yield zero initial gradients and result in no evolution at all.

Figure 4 demonstrates very good overall agreement between the analytical and numerical solutions; the shocks and contact discontinuities are well captured. Although the baryon evolution is decoupled, this test still demonstrates excellent performance of the root finding algorithm.

4.2 Bjorken flow

In this subsection we test BEShydro in Milne coordinates for a transversally homogeneous dissipative fluid undergoing longitudinally boost-invariant Bjorken expansion [47]. Boost-invariant systems are characterized by space-time rapidity independent macroscopic observables and a flow profile that looks static (i.e. $u^{\mu}=(1,0,0,0)$) in Milne coordinates [44]. In spite of experimental evidence for longitudinal density gradients, there are strong phenomenological indications that near mid-rapidity a longitudinally boost-invariant flow profile is a good approximation for relativistic heavy-ion collisions at $\sqrt{s}\gtrsim 100$ GeV/nucleon (see, e.g., [44]). The additional assumption of transverse homogeneity, however, is clearly unrealistic, given the finite transverse size of the colliding nuclei. Still, it provides a useful test bed because the resulting independence of all macroscopic quantities from all three spatial dimensions simplifies the dissipative hydrodynamic evolution equations to a set of coupled ordinary differential equations which can be solved with Mathematica:

[TABLE]

where $\pi\equiv-\tau^{2}\pi^{\eta\eta}$ has been introduced.

In the BEShydro simulation the numbers of cells in $(x,y,\eta_{s})$ directions are all set to 1.111111Of course, one can also set the number of cells in $(x,y)$ larger than 1 but to ensure transverse spatial homogeneity quantities like ${\cal E}$ and ${\cal P}_{0}$ should be the same in all cells. With boost-invariant and transversally homogeneous initial conditions, the numerical results from the (3+1)D BEShydro code can be tested against a solution of the coupled ODEs (106-110) from a separate Mathematica code. As for the Riemann problem, the baryon evolution decouples from the rest of the system if an EoS ${\cal P}_{0}={\cal P}_{0}({\cal E})$ is used; in this case, even with a nonzero baryon diffusion current in longitudinal direction, the Bjorken expansion remains unchanged. Although longitudinal boost-invariance and transverse homogeneity do not allow any gradients of the chemical potential that could drive a baryon diffusion current, we can then still test the baryonic sector of the code by initiating it with nonzero initial values for the net baryon density and diffusion current. These out-of-equilibrium initial values will then relax according to Eqs. (109,110), without affecting the Bjorken flow profile.

By using EOS2 (the Wuppertal-Budapest EoS at $\mu=0$) which features a non-zero interaction measure, we can also test the evolution of the bulk viscous pressure $\Pi$ which is propagated by Eq. (107) as an additional dissipative degree of freedom [37]. Noting that Eq. (107) involves a rather complex parametrization of the transport coefficients (41)-(44), we emphasize that this test is indeed non-trivial.

As Fig. 5 shows, the agreement of BEShydro with the semi-analytic Mathematica solution is excellent. The evolution of the baryon diffusion current follows the exact decay very precisely over six orders of magnitude. This would not be possible if the root finding algorithm did not perform with high accuracy.121212As for the Riemann problem, this test again checks this property only for the case $\partial{\cal P}_{0}/\partial{\cal N}=0$.

Still, because of boost-invariance and transverse homogeneity many terms in the full set of evolution equations vanish in this example. The absence of any kind of transverse expansion in this test is particularly worrisome. This question will be addressed next with the “Gubser test”.

4.3 Gubser flow

Gubser flow [92, 93] describes conformally symmetric systems that, in addition to longitudinally boost-invariant Bjorken flow, undergo at the same time strong azimuthally symmetric (“radial”) transverse flow. Contrary to heavy-ion collisions where transverse flow is initially zero and only generated after the collision in response to transverse pressure gradient, in Gubser flow the transverse flow exists at all times (i.e. even on a hypersurface corresponding to some very early “initial” longitudinal proper time $\tau_{0}$).

4.3.1 Gubser coordinates

Gubser flow originates from an ingenious symmetry that, like Bjorken symmetry, makes the flow appear static in an appropriately chosen set of coordinates called “Gubser coordinates” which map Minkowski space onto a 3-dimensional de Sitter space times a line, dS${}_{3}\otimes\mathbb{R}$ [92, 93]. As a result of this symmetry, macroscopic quantities do not depend on any of the space-like coordinates but only on the time-like coordinate in this system, and the dissipative hydrodynamic equations again reduce to a set of coupled ODEs in that time coordinate, $\rho\in\mathbb{R}$.

To introduce the Gubser coordinates we first rescale the invariant distance measure of Minkowski space in Milne coordinates with a Weyl transformation131313To make our equations readily comparable with those in the by now vast literature on Gubser flow we temporarily (i.e. in this subsection only) switch our metric signature convention to the mostly-plus metric, i.e. $g^{\mu\nu}=\textrm{diag}(-1,1,1,1)$ in Cartesian coordinates.

[TABLE]

Next we perform the coordinate transformation $x^{\mu}=(\tau,r,\phi,\eta_{s})\to\hat{x}^{\mu}=(\rho,\theta,\phi,\eta_{s})$,141414Here $r\equiv\sqrt{x^{2}+y^{2}}$ and $\phi\equiv\tan^{-1}(x/y)$. All quantities expressed as functions of Gubser coordinates are made unitless by scaling them with the appropriate powers of the Milne time $\tau$ and labeled with a hat. by introducing [92, 93]

[TABLE]

where $q$ is an arbitrary energy scale that defines the physical size of the system (the solution is invariant under a common rescaling of $q$, $\tau$ and $r$ such that $qr$ and $q\tau$ remain unchanged). In these coordinates the Weyl-rescaled invariant distance measure becomes

[TABLE]

with the metric

[TABLE]

A system that appears static in the coordinates $\hat{x}^{\mu}=(\rho,\theta,\phi,\eta_{s})$, i.e. has flow velocity $\hat{u}^{\mu}=(1,0,0,0)$, is said to exhibit Gubser flow in Minkowski space.

To map quantities expressed in Gubser coordinates back to Milne coordinates in Minkowski space one uses metric rescaling [93] and the definitions (112,113), for example

[TABLE]

With these transformation rules the Gubser flow profile can be expressed in Milne coordinates through the components

[TABLE]

where $\kappa(\tau,r)$ is the transverse flow rapidity, corresponding to the transverse flow velocity

[TABLE]

Note the transverse flow components are azimuthally symmetric. This flow is dictated by symmetry so it applies to both ideal and dissipative fluids whose thermodynamic functions have Gubser symmetry (i.e. depend only on $\rho$ when expressed in Gubser coordinates). Different initial conditions for the hydrodynamic quantities and different transport coefficients yield different $\rho$ dependencies for their evolution, translating into different characteristics $r(\tau)$ and different $(\tau,r(\tau))$ profiles when expressed in Minkowski space coordinates.

To ensure invariance of the hydrodynamic equations under the Weyl transformation (111) the energy momentum tensor must be traceless. This means that a conformal EoS must be used and the bulk viscosity and bulk viscous pressure must be set identically to zero.

4.3.2 Gubser flow with baryon diffusion

In this work we extend the existing semi-analytical solutions for conformal Israel-Stewart hydrodynamics with Gubser flow [97] to systems with a non-zero baryon diffusion current. One could argue that the longitudinal reflection symmetry under $\eta_{s}\to-\eta_{s}$ in the Gubser symmetry indicates that the baryon diffusion current should be zero [39]. However, when the conformal EoS is used, as required by the Weyl invariance of the hydrodynamics, the baryon evolution decouples from the rest of the system, and the baryon density and diffusion current evolve as background fields. This means that a non-zero baryon diffusion current does not modify the Gubser flow profile, and the numerical results for the evolution of the baryon diffusion current can be tested by comparing them to the semi-analytical solutions of the equations of motion obtained from $\hat{u}^{\mu}=(1,0,0,0)$ under Gubser symmetry.

In this subsection we derive these baryon equations of motion in de Sitter space. We rewrite Eqs. (8, 9) and (18, 19) with the mostly-plus metric tensor and apply the Gubser flow profile $\hat{u}^{\mu}=(1,0,0,0)$ and the de Sitter metric (115) to obtain151515Similar to Bjorken flow, the shear stress for Gubser flow has only one independent component for which we take $\hat{\pi}^{\eta\eta}$. The other non-vanishing components $\hat{\pi}^{\theta\theta}$ and $\hat{\pi}^{\phi\phi}$ are related to $\hat{\pi}^{\eta\eta}$ by tracelessness (which gives $\hat{\pi}^{\eta}_{\eta}=-\hat{\pi}^{\theta}_{\theta}-\hat{\pi}^{\phi}_{\phi}$) and azimuthal symmetry (which implies $\hat{\pi}^{\eta}_{\eta}=-2\hat{\pi}^{\theta}_{\theta}=-2\hat{\pi}^{\phi}_{\phi}$ for evolution with azimuthally symmetric initial conditions).

[TABLE]

(Note that $\hat{\theta}=2\tanh\rho$ is the scalar expansion rate for Gubser flow.) For the transport coefficients we use the same parametrization as described before in Fig. 5 (see also the caption of Fig. 6). The transformation rules for the baryon density and diffusion current are

[TABLE]

Given initial conditions $\hat{\cal E}(\rho_{0})$, $\hat{\pi}^{\eta\eta}(\rho_{0})$, $\hat{\cal N}(\rho_{0})$, and $\hat{n}^{\eta}(\rho_{0})$ at, say, $\rho_{0}=0$, obtaining $\hat{\cal P}_{0}(\rho_{0})$ and $\hat{T}(\rho_{0})$ with the help of EOS1, Eqs. (125-128) can be solved with Mathematica in Gubser coordinates and then translated into Milne coordinates $(\tau,x,y,\eta_{s}=0)$ using the transformation rules (116-119) and (129-130). This semi-analytic solution can then be used to obtain initial conditions for BEShydro on an initial proper time hypersurface $\tau_{0}$ which are then further evolved with the (3+1)-dimensional BEShydro code. Specifically, BEShydro requires initial data at $\tau_{0}$ for

[TABLE]

on the computational $(x,y)$ grid (due to longitudinal boost-invariance these are only required at $\eta_{s}=0$). All remaining hydrodynamic components are either zero or can be obtained from the above by symmetry. For example, all shear stress components $\pi^{\mu\nu}(\tau_{0},x,y,\eta_{s}=0)$ for BEShydro can be obtained from the semi-analytic solution $\hat{\pi}^{\eta\eta}(\rho)$ by using tracelessness and azimuthal symmetry,

[TABLE]

(with all other Gubser components being zero by symmetry), followed by

[TABLE]

This gives, for example, at $(\tau{=}\tau_{0},\eta_{s}{=}0)$

[TABLE]

where the value of $\rho$ depends on the transverse grid point $(x,y)$.161616Obviously, whenever the transport coefficients are changed in BEShydro, the semi-analytic solution must be recomputed accordingly for comparison, also because the full exact solution (not just its initial conditions in Gubser coordinates) is required to obtain initial conditions for the BEShydro code in Milne coordinates.

In Fig. 6 we compare the BEShydro output with the semi-analytic Gubser solution for the default setup described in the figure caption. One observes excellent agreement. Owing to the non-trivial transverse expansion of Gubser flow it allows to test additional source terms in the BEShydro evolution equations when compared with the Bjorken flow test. We can also use it to study the performance of the root-finding algorithm in BEShydro.

As discussed in Sec. 3.3, two different methods are used in complementary ranges of the flow velocity separated by the critical value $v=0.563624$ or, equivalently, $u^{\tau}=1.21061$. Fig. 6(e) shows that the root-finding algorithm works equally well on both sides of the critical value of $u^{\tau}$. We also checked the precision and relative speed of convergence of the Newton-Raphson and modified iteration schemes described in Sec. 3.3. For both methods, the maximum number of iterations is set to 100 (which is never reached), and the root finding stops when the relative error $|v_{i+1}-v_{i}|/v_{i+1}<10^{-6}$ or $|u^{\tau}_{i+1}-u^{\tau}_{i}|/u^{\tau}_{i+1}<10^{-4}$. For the conformal EoS used here we found both methods to converge about equally well to the same result (within the specified uncertainty), with the modified iteration scheme being about 15% faster than the Newton-Raphson method.171717During the early evolution stages the Newton-Raphson method converges somewhat faster but at later times the modified interaction scheme is found to be more efficient.

4.4 Comparison to other codes

Owing to their shared assumption of longitudinal boost-invariance, none of the tests described in the preceding subsections addresses the performance of the BEShydro code in describing the expansion along the rapidity direction. To remedy this we have compared BEShydro output for a transversally homogeneous system undergoing arbitrary longitudinal expansion without transverse expansion with the results from an independent (1+1)D hydrodynamic code developed by Monnai [94]. To be able to study baryon number transport, EOS4 is used in both codes. A similar comparison was also made in Ref. [39] to test the performance of MUSIC. Rather than directly using Monnai’s code [94] we compare our BEShydro results with those reported in the comparison [39] with MUSIC. In this sense, the following test is also a code comparison with MUSIC. Qualitatively consistent results from an earlier BEShydro study were already reported in [27].

We here focus on the baryon number evolution in the longitudinal direction, with non-vanishing longitudinal gradients of $\mu/T$, by setting bulk and shear stresses to zero.181818In principle, since shear stress is known to affect the evolution of the energy density it might be expected to also change the baryon number flow; this is an interesting physics question which we postpone for a separate study. Specifically, we check how BEShydro handles the evolution equation (30) for baryon number diffusion,

[TABLE]

together with the net baryon conservation law (28). For the baryon transport coefficients we choose $\kappa_{n}=0.2\,{\cal N}/\mu$ and $\tau_{n}=0.2/T$.

In Fig. 7 we show a comparison of the distributions in space-time rapidity $\eta_{s}$ of energy density ${\cal E}$ and net baryon density ${\cal N}$ at four different times. Broken (continuous) lines show the results from BEShydro (MUSIC [39]); the agreement between these two codes is very good. In addition to testing the longitudinal dynamics this comparison also shows that the root-finding algorithm works correctly with a realistic EoS  ${\cal P}_{0}({\cal E},{\cal N})$ that depends on both energy and net baryon density.

4.5 Tests summary

We briefly summarize which parts of the BEShydro code were tested with the test protocols described in this section. As described in Sec. 3, the same RK-KT algorithm is applied for solving the equations of motion for all hydrodynamical variables propagating in the system, and the same root-finding algorithm is used for all equations of state, whether they depend on baryon density or not. The observed good performance in propagating all hydrodynamical variables indicates the efficiency of the RK-KT algorithm, and the proper evolution of baryon density and baryon diffusion validates the root-finding method in systems with non-zero baryon currents.

The Riemann problem in ideal hydrodynamics with EOS1 (Sec. 4.1) shows the ability of the RK-KT algorithm in BEShydro to capture shocks and contact discontinuities. The (0+1)D Bjorken expansion with EOS2 (Sec. 4.2) tests the programming of the equations of motion, especially for the shear and bulk components, including the non-trivial parametrization of the transport coefficients, and the root-finding algorithm with non-zero but decoupled baryon density and diffusion current. The Gubser flow test with EOS1 (Sec. 4.3) provides extra validation in situations with strong transverse expansion featuring large temporal and transverse gradients. The comparison with MUSIC and Monnai’s (1+1)D code (Sec. 4.4) at finite baryon density with the realistic EOS4 validates the longitudinal dynamics of density and baryon diffusion without the simplification of longitudinal boost-invariance, as well as the root-finding algorithm with non-zero baryon density and baryon diffusion, including the bilinear interpolation of the EoS tables. The figures shown in this section demonstrate that BEShydro passes all these tests without struggle.

5 Baryon diffusion in an expanding QGP

In this section we illustrate the evolution of energy and baryon number in an expanding QGP with realistic “bumpy” initial conditions, by visualizing the evolution of the corresponding densities in the transverse plane at $\eta_{s}=0$.191919For a similar earlier study with smooth, ensemble-averaged initial conditions see [39]. This generalizes many similar visualizations made in the past for systems without conserved charges. Since the physics of initial-state fluctuations along the longitudinal direction is still not very well explored (for a few examples see, Refs. [52, 94, 39]), we here use smooth longitudinal initial conditions and refrain from showing the (mostly uninteresting) evolution along the beam direction.

Following Refs. [98, 27, 39] we use a 3-dimensional initial condition at non-zero baryon density which extends a transverse profile obtained from the Monte Carlo (MC) Glauber model [99] into the longitudinal direction with the following prescription:

[TABLE]

where $T_{A/B}(x,y)$ are the transverse profiles of the right- and left-moving nuclei from the MC-Glauber model and ${\cal E}_{A/B}(\eta_{s})$, ${\cal N}_{A/B}(\eta_{s})$ are the corresponding longitudinal profiles for the energy and net baryon density, respectively [98, 27, 39]. ${\cal E}_{0}$ is a normalization factor which can be tuned to reproduce the final multiplicity while ${\cal N}(\tau_{0},x,y,\eta_{s})$ is normalized to the total number of participant baryons [39].

In Fig. 8, we show the transverse distributions at $\eta_{s}=0$ of the energy density (upper panels) and net baryon density (lower panels), at the hydrodynamic starting time $\tau_{0}=0.5$ fm/$c$ (left column) and later at $\tau=5$ fm/$c$, for evolution with only shear stress turned on (middle column) and both shear stress and baryon diffusion turned on (right column). By comparing the middle and right columns we note that baryon diffusion leaves no pronounced signatures in the evolution of the energy density but smoothes out gradients in baryon density. The authors of Refs. [98, 27, 39] came to a similar conclusion for the evolution in the longitudinal direction.

6 Conclusions and outlook

In this work we have described the structure and performance of a new code called BEShydro describing the (3+1)-dimensional space-time evolution of hot and dense matter created in relativistic heavy-ion collisions using second-order dissipative relativistic fluid dynamics. BEShydro differs from most other publicly available algorithms by evolving, together with energy and momentum, a conserved current related to net baryon number, including its dissipative diffusion current, and it evolves the medium with an equation of state that depends on both the energy and net baryon densities. The generalization of the BEShydro algorithm to the simultaneous propagation of multiple conserved charge currents [41] should be a straightforward task for the near future. A dynamical initialization module for BEShydro that accounts for the non-zero temporal extension of the energy deposition process in heavy-ion collisions at collision energies probed in the Beam Energy Scan (BES) campaign at the Relativistic Heavy-Ion Collider, along the lines discussed in Refs. [26, 27], is under construction and will soon be added to the BEShydro package. The code can also be plugged in the JETSCAPE framework [100] as a hydrodynamic evolution module.

After briefly describing the physics to be addressed with BEShydro simulations, we presented in detail the equations BEShydro solves, the transport coefficients and equations of state it uses, and the numerical schemes it employs, including the root-finding algorithm for finding the flow velocity field at each time step and the regulation schemes used to regulate overly large dissipative flows caused by numerical or anomalously large physical fluctuations. The performance of the code was verified with high precision using a series of standard tests involving comparison with analytically or semi-analytically known solutions for problems of reduced dimensionality, characterized by additional symmetries that are usually not respected in real life situations but can be exploited for code verification. The code is distributed together with this suite of verification test protocols, thus enabling the user to check the continued accuracy of the code after changing or generalizing some of its parts. Finally, we presented a simple example illustrating the practical effects of baryon number diffusion on the evolution of energy and net baryon density for a collision between copper (Cu) nuclei, demonstrating the smoothing effects of baryon diffusion on large gradients of the net baryon density in the transverse plane.

Acknowledgements

The authors thank Dennis Bazow, Derek Everett, Mike McNelis, Long-Gang Pang, Dirk Rischke, Chun Shen, Leonardo Tinti and Gojko Vujanovic for fruitful discussions. This work was supported in part by the U.S. Department of Energy (DOE), Office of Science, Office for Nuclear Physics under Award No. DE-SC0004286 and within the framework of the BEST Collaboration, and in part by the National Science Foundation (NSF) within the framework of the JETSCAPE Collaboration under Award No. ACI-1550223. LD appreciates the kind hospitality of the European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*) and the MIT Center for Theoretical Physics during the completion of this work. UH’s stay at the Institut für Theoretische Physik of the J. W. Goethe-Universität was supported by the Alexander von Humboldt Foundation through a Research Prize. Computing resources were generously provided by the Ohio Supercomputer Center [101] and by the Open Science Grid [102, 103], which is supported by the National Science Foundation, award 1148698, and the U.S. Department of Energy’s Office of Science.

Appendix A Conservative form of the evolution equations

In this Appendix, we recast Eqs. (24-28) and (29-31) in the conservative form (66). The procedure follows Refs. [37, 15], adding here the equations for baryon evolution. The equations reproduced here are written in a form that facilitates direct comparison with the BEShydro code.

We start with the conservation laws in Eqs. (24-28). With the scaled flow velocities $v^{i}=u^{i}/u^{\tau}$ $(i=x,y,\eta_{s})$, we can write down the following constituent relations for the components of $T^{\mu\nu}$ and $N^{\mu}$:

[TABLE]

Inserting these into the conservation laws one obtains

[TABLE]

with the following source terms for $T^{\tau\tau}$:

[TABLE]

for $T^{\tau x}$:

[TABLE]

for $T^{\tau y}$:

[TABLE]

for $T^{\tau\eta}$:

[TABLE]

and for $N^{\tau}$:

[TABLE]

Considering $d\equiv u^{\mu}\partial_{\mu}$, the relaxation equations (29-31) for the dissipative flows can be written as

[TABLE]

where we used $\partial_{i}v^{i}\equiv\partial_{x}v^{x}+\partial_{y}v^{y}+\partial_{\eta}v^{\eta}$. The source terms in Eqs. (167-169) are given by

[TABLE]

Here the $I$-terms and $G$-terms are defined in Eqs. (29)-(31).

Appendix B Explicit form

In this Appendix we provide explicit forms of some expression occurring in the baryon evolution equations. In Milne coordinates the terms $G^{\mu}_{n}=u^{\alpha}\Gamma^{\mu}_{\alpha\beta}n^{\beta}$ in Eq. (30) evaluate to

[TABLE]

To calculate the LRF gradient of $\mu/T$ in Eq. (30) numerically we use

[TABLE]

and work out the partial derivatives $\partial_{\mu}(\mu/T)$ in the computational frame numerically from the EoS tables. Finally, the last two terms in Eq. (35) can be expressed as

[TABLE]