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





a]Lipei Duauthor
a,b]Ulrich Heinzemail
[author] Email: du.458@osu.edu (corresponding author)
[email] Email: heinz.9@osu.edu
[a]Department of Physics, The Ohio State University, Columbus, OH 43210-1117, USA
[b]Institut für Theoretische Physik, J. W. Goethe Universität, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany



    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  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.


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







PROGRAM SUMMARY




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

Authors: Lipei Du, Ulrich Heinz 

Program Title: 

Licensing provisions: GPLv3 

Programming language: C++  

Computer: Laptop, desktop, cluster  

Operating system: GNU/Linux distributions, Mac OS X 

Memory usage: For a 121×121×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, 

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×121 uses 0.32 sec/time step, and on a 3-dimensional grid 121×121×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.








§ INTRODUCTION



 Over the last decade, second-order dissipative relativistic fluid dynamics (RFD) <cit.> 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 <cit.>. 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 <cit.> 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 <cit.>, and it must be coupled to a microscopic kinetic description of the hadronic rescattering stage at the end of the fireball evolution <cit.>, 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 <cit.>, 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 <cit.> and FAIR <cit.>.

We here describe a new (3+1)-dimensional dissipative RFD code, which we call , that includes the evolution of the net baryon charge and diffusion currents.[
    The code is open source and can be freely downloaded from <https://github.com/LipeiDu/BEShydro>]

 evolved out of the CPU version of the code GPU-VH <cit.>, by adding a number of additional dissipative terms <cit.> 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  is the latest version of MUSIC <cit.>, while other codes (e.g. <cit.>) so far ignore the evolution of the net baryon diffusion current or, in some cases, even that of the net baryon charge.  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 .
We start from the most general form of the equations solved by  in Sec. <ref>, describe the numerical scheme used by the code in  Sec. <ref>, and include a number of validation tests in Sec. <ref>, 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  evolves only a single conserved charge (net baryon number); future generalization to include strangeness and isospin <cit.> is expected to be conceptually straightforward. 





§ RELATIVISTIC HYDRODYNAMICS




In this section the physics ingredients of  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.




 §.§ Equations of motion






  §.§.§ Conservation laws



Hydrodynamics is a macroscopic theory describing the space-time evolution of the 14 components of the energy-momentum tensor T^μν(x) and the net (baryon) charge current N^μ(x).[
     For 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 <cit.>:


    d_μT^μν   ≡   1/√(g)∂_μ(√(g)T^μν)+Γ^ν_μλT^μλ=0 ,
    
    
    d_μN^μ   ≡   1/√(g)∂_μ(√(g)N^μ)=0 .


Here d_μ (μ=0,1,2,3) stands for the covariant derivative in a general system of space-time coordinates, with metric tensor g^μν defined with negative signature (“mostly minus” convention (+,-,-,-)), g≡-(g_μν), and the Christoffel symbols (see, e.g., <cit.>)


    Γ^μ_αβ≡1/2 g^μν(∂_βg_αν
             +∂_αg_νβ
             -∂_νg_αβ)
        =Γ^μ_βα .



The 14 independent components of T^μν and N^μ are more physically defined in terms of the hydrodynamic decomposition of these tensors <cit.>,


    T^μν   =    u^μu^ν-(+Π)Δ^μν+π^μν ,
    
        N^μ   =    u^μ+^μ .


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


    T^μνu_ν= u^μ ,


and specifies the local rest frame (LRF) of the fluid at point x (the so-called “Landau frame”). The tensors u^μu^ν and Δ^μν≡ g^μν - u^μu^ν then are projectors on the temporal and spatial directions in the LRF.  and  are the energy and net baryon density in the LRF which can be obtained as the following projections of T^μν and N^μ:


    = u_μ T^μν u_ν ,     = u_μ N^μ .


From these, the local equilibrium pressure  is obtained through the EoS =(,). The shear stress π^μν, the bulk viscous pressure Π, and the baryon diffusion current n^μ are dissipative flows describing deviations from local equilibrium.

Using the decomposition (<ref>,<ref>) the conservation laws (<ref>,<ref>) can be brought into the physically intuitive form <cit.> 


    D   =    -θ -∇_μ n^μ ,
    
    
     D   =    -(++Π)θ + π_μνσ^μν ,
    
    
      (++Π)  Du^μ   =   ∇^μ(+Π)
        - Δ^μν∇^σπ_νσ + π^μν Du_ν .


Here D=u_μ d^μ denotes the time derivative in the LRF, θ=d_μ u^μ is the scalar expansion rate, ∇^μ=∂^⟨μ⟩ (where generally A^⟨μ⟩≡Δ^μνA_ν) denotes the spatial gradient in the LRF, and σ^μν = ∇^⟨μu^ν⟩ (where generally B^⟨μν⟩≡Δ^μν_αβ B^αβ, with the traceless spatial projector Δ^μν_αβ≡1/2 (Δ^μ_αΔ^ν_β + Δ^ν_αΔ^μ_β) - 1/3Δ^μνΔ_αβ) is the shear flow tensor. While these equations clearly exhibit the physics in the LRF (which varies from point to point),  solves the conservation laws (<ref>,<ref>) in a fixed global computational frame. Their explicit form in the global frame is discussed in Sec. <ref>.

The 5 conservation laws (<ref>-<ref>) are sufficient to determine the energy and baryon density,  and , together with the 3 independent components of the flow velocity u^μ, as long as the shear stress π^μν, the bulk viscous pressure Π, and the baryon diffusion current n^μ vanish.[
     Note that the shear stress is traceless, π^μ_μ=0, and both π^μν and n^μ have only spatial components in the LRF, u_μπ^μν=π^μν u_ν = u_μ n^μ =0. π^μν, Π, and n^μ 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 <cit.> 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 μ, the equilibrium pressure (,)=(T,μ) is recognized as the grand-canonical thermodynamic potential for a system with temperature T(,) and chemical potential μ(,). In principle, T and μ 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  the driving force for net baryon number diffusion is formulated in terms of the gradient of μ/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 μ since they are defined as response functions of the thermal equilibrium system described by the potential (T,μ). Different versions of the EoS (T,μ) or (,) used in  will be described in Sec. <ref>.




  §.§.§ Evolution equations for the dissipative flows




In  the dissipative flows are evolved with DNMR theory <cit.>. While the equations of motion in this theory are derived from the Boltzmann equation which is applicable only to weakly coupled systems <cit.>, 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 <cit.>. 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:


    τ_ΠDΠ+Π   =   Π_NS +𝒥+𝒦+ℛ ,
    τ_n(Dn)^⟨μ⟩+n^μ   =    n^μ_NS+𝒥^μ+𝒦^μ+ℛ^μ ,
    τ _π(Dπ)^⟨μν⟩+π ^μν   =   π^μν_NS+𝒥^μν+𝒦^μν+ℛ^μν .


Here (Dn)^⟨μ⟩≡Δ^μν Dn_ν and (Dπ)^⟨μν⟩≡Δ^μν_αβ Dπ^αβ, ensuring that all terms are purely spatial in the LRF and, where applicable, traceless. τ _Π, τ_n, and τ _π are the relaxation times for Π, ^μ, and π^μν, respectively. They control how fast the dissipative flows relax to their Navier-Stokes limits <cit.>:


    Π_NS   =    -ζθ ,
    
        n^μ_NS   =   κ_n ∇^μ(μ/T) ,
    π ^μν_NS   =    2ησ ^μν ,


where ζ, κ_n, and η 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 θ, the spatial gradient of μ/T in the LRF, ∇^μ(μ/T), and the shear flow tensor  σ ^μν, respectively, which drive the system away from local equilibrium. 

The scalar, vector and tensor source terms on the r.h.s. of Eqs. (<ref>-<ref>), 𝒥, 𝒦, ℛ, 𝒥^μ, 𝒦^μ, ℛ^μ, 𝒥^μν, 𝒦^μν and ℛ^μν 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 <cit.>, the Navier-Stokes and 𝒦 terms on the r.h.s. of Eqs. (<ref>-<ref>) are of first and second order in the Knudsen number(s), respectively, the  𝒥 terms are of order Knudsen number times inverse Reynolds number, and the ℛ terms are of second order in the inverse Reynolds number(s). Their explicit expressions can be found in Ref. <cit.>. Following the arguments in Ref. <cit.>, we here include only a subset of the 𝒥 terms. As  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


    τ_ΠDΠ+Π   =   
            -ζθ -δ_ΠΠΠθ
            +λ_Πππ^μνσ_μν ,
    
    τ _n D n^μ+n^μ   =   κ_n ∇^μ(μ/T) -τ_n n_νω ^νμ-δ_nnn^μθ
       -   λ _nnn_νσ ^μν - τ_n n^ν u^μ D u_ν ,
    
    τ _πDπ^μν+π ^μν   
            =   2ησ ^μν
            +2τ_ππ^⟨μ_λω^ν⟩λ
            -δ_πππ^μνθ
       -   τ_πππ^λ⟨μσ^ν⟩_λ
            +λ_πΠΠσ^μν
       -   τ_π( π ^λμu^ν+π ^λνu^μ) Du_λ ;


here ω^μν=1/2(∇^μ u^ν-∇^ν u^μ) is the vorticity tensor. The additional transport coefficients δ_ΠΠ, λ_Ππ, τ_n etc. will be discussed in Sec. <ref>. In Eqs. (<ref>) and (<ref>) we removed the transversality constraints on the l.h.s. by using footnote <ref> and


    (Dn)^⟨μ⟩ = Δ^μ_ν D n^ν = D n^μ + u^μ n^ν D u_ν ,


as well as its analog for π^μν, moving the extra terms as additional source terms to the r.h.s.

It is worth pointing out that in Eqs. (<ref>)-(<ref>) we have followed Ref. <cit.> in ignoring terms describing the direct influence of baryon diffusion, n^μ, on the evolution of the shear and bulk viscous stresses, π^μν and Π. Baryon evolution still affects the evolution of the system indirectly through the EoS. In this approach it has been shown <cit.> 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.




  §.§.§ 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^μ=(τ,x,y,η_s) which are best adapted to the relativistic collision kinematics and subsequent almost boost-invariant longitudinal flow pattern <cit.>. In terms of Cartesian coordinates (t,x,y,z) the longitudinal proper time τ and space-time rapidity η_s are defined as


    τ=√(t^2-z^2) ,    η_s=1/2ln(t+z/t-z) .


The mid-rapidity point z=η_s=0 at τ=0 defines the collision point in the global (computational) frame. In Milne coordinates the metric tensor is


    g^μν =diag(1,-1,-1,-1/τ^2) ,


the fluid four-velocity is u^μ =(u^τ,u^x,u^y,u^η), and the four-derivative is ∂_μ=(∂_τ,∂_x,∂_y,∂_η).[
        In all sub- and superscripts η is short for η_s.]

The metric has the following non-vanishing Christoffel symbols:


    Γ^η_τη=Γ^η_ητ=1/τ ,    Γ^τ_ηη=τ .


Plugging them into Eqs. (<ref>,<ref>) we obtain the conservation laws in Milne coordinates:


    ∂_μT^μτ    =   -1/τ(T^ττ+τ^2T^ηη)  , 
    ∂_μT^μ x    =   -1/τT^τ x  , 
    
    ∂_μT^μ y   =   -1/τT^τ y  , 
    
    ∂_μT^μη   =   -3/τT^τη ,
    
    ∂_μ N^μ   =   -1/τ N^τ .


Introducing the convective time derivative d ≡ u^μ∂_μ, the relaxation equations can be written as


    dΠ   =   -ζ/τ_Πθ-Π/τ_Π-I_Π ,
    
        d ^μ   =   κ_n/τ_∇^μ(μ/T) - ^μ/τ_ - I_^μ 
        - G_^μ ,
        
    
        dπ^μν   =   2η/τ_πσ^μν-π^μν/τ_π-I^μν_π-G^μν_π ,


with the shorthand notations G^μν_π≡ u^αΓ^μ_αβπ^βν + u^αΓ^ν_αβπ^βμ and G^μ_n = u^αΓ^μ_αβ n^β for the geometrical source terms obtained when splitting the covariant LRF time derivative D in Eqs. (<ref>)-(<ref>) into the convective time derivatives d and a remainder (for example, D n^μ≡ u^α d_α n^μ = u^α (∂_α n^μ + Γ^μ_αβ n^β)). The I-terms are explicitly


    I_Π   ≡   δ_ΠΠ/τ_ΠΠθ-λ_Ππ/τ_Ππ^μνσ_μν ,
    
    I^μ_n    ≡    I_1^μ+ δ _nn/τ _n I_2^μ+ I_3^μ+ λ _nn/τ _n I_4^μ , 
    
    I^μν_π   ≡    I^μν_1
    +δ_ππ/τ_πI^μν_2
    -I^μν_3
    +τ_ππ/τ_πI^μν_4
    -λ_πΠ/τ_πΠσ^μν ,


with 


    I_1^μ = u^μ n^ν Du_ν,   
        I_2^μ = n^μθ,  
        I_3^μ = n_νω^νμ,  
        I_4^μ = n_νσ^νμ;
    
          
        I_1^μν = (u^μπ^νλ+u^νπ^μλ) Du_λ,   
        I_2^μν = θπ^μν,   
        I_3^μν = ω^μ_ λπ^λν
                        +ω^ν_ λπ^λμ,    
          
        I_4^μν = 1/2(π^μλσ^ ν_λ
             +π^νλσ^ μ_λ) 
       -1/3Δ^μνπ^αβσ_βα.



The conservation laws (<ref>-<ref>) together with the dissipative transport equations (<ref>-<ref>) constitute the equations of motion of the relativistic hydrodynamic system encoded in . Next we will discuss the transport coefficients appearing in these equations.




 §.§ 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. <ref>), 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,  implements the transport coefficients from Ref. <cit.> which starts from the Boltzmann equation in Relaxation Time Approximation (RTA) and employs the 14-moment approximation <cit.> for a one-component gas of Boltzmann particles with non-zero but small mass m≪ T. In all expressions we only keep the leading non-zero term in powers of m/T≪ 1. This is motivated by the approximate masslessness of the microscopic quark-gluon degrees of freedom that make up the fluid described by .




  §.§.§ Shear stress tensor



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


    η/τ_π   =   +/5 , 
    
    δ_ππ/τ_π = 4/3 ,     τ_ππ/τ_π   =   10/7 ,     λ_πΠ/τ_π = 6/5 .


Following <cit.> we express shear viscous effects in terms of the kinematic shear viscosity 


    η̅= η T/+ .


For μ=0 the kinematic shear viscosity reduces to the specific shear viscosity η/, where  denotes the entropy density, but it differs greatly from η/ at large net baryon densities, and this can lead to significant differences in the hydrodynamic flow patterns <cit.>. Parametrizations of η̅ as a function of T and μ are discussed in <cit.>; in  the default setting for η̅ is a constant, η̅=0.2. Given η̅, Eq. (<ref>) is used to calculate τ_π=5η̅/T, and Eqs. (<ref>) to obtain the remaining transport coefficients. 




  §.§.§ Bulk viscous pressure



To evolve the bulk viscous pressure we take the transport coefficients 


    ζ/τ_Π   =   15(1/3-c^2_s)^2(+) , 
    
    δ_ΠΠ/τ_Π   =   2/3 ,    λ_Ππ/τ_Π=8/5(1/3-c^2_s) ,


where c_s is the speed of sound in the medium (see Eq. (<ref>)). Similar to the shear viscosity, we start from a parametrization as a function of T and μ of the kinematic bulk viscosity


    ζ̅= ζ T/+


and determine from it the bulk relaxation time τ_Π using Eq. (<ref>) and the remaining transport coefficients using Eqs. (<ref>). 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 <cit.>:


    ζ̅ =
      
      A_0 + A_1 x + A_2 x^2 ,    0.995 T_c ≥ T ≥ 1.05 T_c , 
      
    λ_1 exp [-(x-1)/σ_1] + λ_2 exp [-(x-1)/σ_2]+0.001 ,
            T > 1.05 T_c , 
      
    λ_3 exp [(x-1)/σ_3] + λ_4 exp [ (x-1)/σ_4]+0.03 ,
            T < 0.995 T_c ,


with x = T/T_c and fitted parameters 


    A_0=13.45 ,   A_1=27.55 ,   A_2=-13.77 , 
          λ_1=0.9 ,  λ_2=0.25 ,  λ_3=0.9 ,  λ_4=0.22 ,
          σ_1=0.025 ,  σ_2=0.13 ,   σ_3=0.0025 ,  σ_4=0.022 .



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 <cit.>. Bulk viscous effects are expected to be strongest near the quark-hadron phase transition <cit.>, especially near the QCD critical point where critically enhanced contributions associated with critical slowing-down play a key dynamical role
<cit.>.




  §.§.§ 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. <cit.>, we use the coefficients obtained from the Boltzmann equation for an almost massless (m/T≪1) classical gas (for a calculation of the baryon diffusion coefficient κ_n for a massive gas of hadrons see Ref. <cit.>):


    κ_n    =   τ_n[1/3(μ/T)- T/+] ,
    
    δ_nn   =   τ_n ,    λ_nn = 3/5τ_n .


Here τ_n is the relaxation time of the baryon diffusion current in Eq. (<ref>) and parametrized as 


    τ_n = C_B/T,


with a free parameter C_B.[
        Note 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 (<ref>) for the diffusion coefficient κ_n was derived in first-order Chapman-Enskog approximation <cit.> whereas the second-order transport coefficients (<ref>) were obtained in the 14-moment approximation <cit.>. In the limit of small net baryon density (μ→0) the diffusion coefficient κ_n reduces to


    κ_n/τ_n =  T/3μ .



In RTA, where the collision term in the Boltzmann equation is parametrized with a single relaxation time , the relaxation times for Π, π^μν and ^μ are all the same, i.e. τ_Π=τ_π = τ_n =. Here we allow them to be different. If one does, however, impose the constraint τ_n=τ_π=, the baryon diffusion and shear viscosity coefficients can be related as follows <cit.>:


    κ_n T/η = C(μ/T) (π T/μ)^2
        ( T/+)^2 .


In contrast to Eqs. (<ref>,<ref>), this expression takes into account quantum statistics. The function C(μ/T) exhibits a weak dependence on μ/T, interpolating between 5/3 at large μ/T and a somewhat smaller value at small μ/T whose precise magnitude depends on the number of massless degrees of freedom in the gas <cit.>. Note that, since for small μ/T the net baryon density ∝μ, both (<ref>) and (<ref>) yield nonzero baryon diffusion coefficients at zero net baryon density. At large μ/T, the authors of Ref. <cit.> have shown that the ratio κ_n T/η 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, <cit.>). 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 <cit.>. Studying and comparing the baryon number evolution with baryon diffusion coefficients corresponding to a weakly <cit.> or a strongly coupled QGP <cit.> can be interesting <cit.>.




 §.§ Equation of State (EoS)



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


    = (, ) = (T,μ) .


In practice, for the calculation of the transport coefficients and chemical forces we also need the equivalent relations
T(,) and μ(,). We use the term EoS generically for any one of these relations.




  §.§.§ 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 <cit.>. 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., <cit.>). The method works well at temperatures above ∼ 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^μν, - 3 (usually referred to as the “trace anomaly” or “interaction measure”), describing deviations from the conformal EoS. Defining the rescaled dimensionless trace anomaly

    I(T, μ) ≡(T, μ) - 3(T, μ)/T^4 ,

the thermal pressure at zero chemical potential can be written as

    ^LAT(T, 0)/T^4 = ∫^T_0 dT'I^LAT(T', 0)/T' .



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,μ) suffers from a “sign problem” <cit.>, precluding its direct computation with standard Monte-Carlo methods. This problem can be partially circumvented by using standard LQCD techniques to also compute the μ-derivatives of the pressure (T,μ) at μ=0 and construct (T,μ) at non-zero μ from its Taylor series around μ=0:


    ^LAT(T, μ)/T^4 = ^LAT(T, 0)/T^4 + ∑_n=1^n_max c_2n(T)(μ/T)^2n .


The expansion coefficients 


    c_n(T) = .1/n!∂^n (/T^4)/∂(μ/T)^n|_μ=0


are known as the “baryon number susceptibilities” <cit.>.
The computational effort of computing them increases rapidly with their order n; at this time, the Taylor expansion (<ref>) includes terms up to order n_max=3 and, near T=T_c, converges well up to about baryon chemical potentials μ/T ≲ 2 <cit.>. 

At low temperatures T≪ T_c≃ 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 (-3)/T^4 is given as  <cit.>

    I^HRG(T, 0)=∑_m_i≤ m_maxg_i/2π^2∑^∞_k=1(-η_i)^k+1/k(m_i/T)^3K_1(km_i/T),

where particle species with spin-isospin degeneracy g_i and mass m_i smaller than some cut-off m_max can be included, and η_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 <cit.> and SMASH <cit.>); 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. <cit.> the pressure is interpolated as follows:


    (T, μ)/T^4   =   1/2[1-tanh(T-T'(μ)/Δ T')]
        ^HRG(T, μ)/T^4
       +   1/2[1+tanh(T-T'(μ)/Δ T')]
        ^LAT(T, μ)/T^4 .

Here ^HRG and ^LAT are the equilibrium pressures for the hadron resonance gas and from lattice QCD, respectively, T'(μ) is the “switching temperature” and Δ T' controls the width of the “overlap region”. The authors of Ref. <cit.>, on the other hand, interpolate the interaction measure at μ=0, I(T,μ=0) smoothly between T_1=T_c(μ=0)=155 MeV and T_2=180 MeV, using a polynomial interpolation function. 

Once the pressure (T, μ) is given, other thermodynamic quantities can be calculated from thermodynamic identities:

    (T, μ)/T^3   =   1/T^3[∂(T, μ)/∂ T]_μ ,
    
    (T, μ)/T^3   =   1/T^3[∂(T, μ)/∂μ]_T ,
    
    (T, μ)/T^4   =   (T, μ)/T^3-(T, μ)/T^4+μ/T(T, μ)/T^3 ,
    
    
        c^2_s(T, μ)    =   [∂(, )/∂]_ + /+[∂(, )/∂]_ .

The last equation requires expressing  and  through T and μ after taking the derivatives. In practice, the functions (T,μ) and (T,μ) are numerically inverted, and the quantities T, μ, ,, and c^2_s, as well as the two derivatives on the r.h.s. of Eq. (<ref>), are stored in a table on a (,) grid which is interpolated by the hydrodynamic code as needed. 




  §.§.§ Equations of state implemented in 



In  four different equations of state, EOS1 to EOS4, are implemented, for different purposes: at zero chemical potential, we include a massless (conformal) EoS (=3, EOS1) as well as an interpolated LQCD-HRG EoS from the Wuppertal-Budapest collaboration <cit.> (EOS2); at non-zero chemical potential, an appropriately generalized conformal EoS (EOS3) and an interpolated LQCD-HRG EoS from Ref. <cit.> (EOS4) are used.

EOS1 assumes an ideal gas of massless quarks and gluons:


    = 3 = 3[2(N_c^2-1)+7/2N_cN_f]
                          π^2/90T^4 ,


where N_c = 3 and N_f = 2.5 are the numbers of colors and (approximately) massless quark flavors, respectively.[
    Strange 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 μ = 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  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 <cit.> can be used for realistic simulations at vanishing μ, i.e. for heavy-ion collisions at ultra-relativistic collision energies with √(s_NN)≫ 100 GeV, especially near mid-rapidity. More details about EOS2 can be found in Refs. <cit.>.

EOS3 is the generalization of EOS1 to non-zero μ. Starting from the ideal massless parton gas expression <cit.>


    /T^4 = π^2/90[2(N_c^2-1)+∑_f 4N_c(7/8+15/4(μ_f/π T)^2+15/8(μ_f/π T)^4)],







where N_c=3 and the sum goes over massless quark flavors, we simplify it by setting μ_f = μ/3 for all flavors (which is appropriate if only baryon number is considered as a conserved charge): 


    (T,μ)/T^4   =   
      p_0+N_f[1/18(μ/T)^2
         +1/324π^2(μ/T)^4]
         =(T,μ)/3T^4 ,
    
    (T,μ)/T^3   =   
       N_f[1/9(μ/T)
               +1/81π^2(μ/T)^3] ,


with p_0=(16+10.5N_f)π^2/90. We again count strange quarks with a factor 1/2, i.e. we set N_f=2.5 so that for μ=0 (<ref>) reduces to (<ref>). Inverting these functions numerically one obtains the EoS tables used in the hydrodynamic code (see Fig. <ref>).

EOS4 from Ref. <cit.>, 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 μ/T using techniques discussed in Sec. <ref>, allows to study the evolution of systems with non-zero net baryon density. It is plotted in Fig. <ref>. In the code, tabulated values for ,  T and μ/T as functions of (, ) are included. If the code requires the EoS at (, ), , T and μ/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 μ. 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 <cit.> and could be imported into  for future dynamical simulations aiming at helping to locate the QCD critical point.






§ NUMERICAL SCHEME




We now describe the numerical scheme used in  to solve the coupled set of evolution equations (the conservation laws (<ref>-<ref>) together with the dissipative relaxation equations (<ref>-<ref>) and the EoS (<ref>)) discussed in Section <ref>. Initial values for all components of the baryon charge current and energy-momentum tensor are set on a surface of constant longitudinal proper time τ.[
        A 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 <cit.>
        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. <cit.> for additional technical details.




 §.§ The Kurganov-Tadmor algorithm




Using the definition of d ≡ u^μ∂_μ, the equations of motion in Eqs. (<ref>-<ref>), (<ref>-<ref>) can all be written <cit.> in the same first-order flux-conserving form


    ∂_τ q + ∂_x (v^x q) + ∂_y (v^y q) + ∂_η (v^η q) = S_q ,


where v^i ≡ u^i/u^τ (i = x, y, η_s) is the 3-velocity of the fluid, the conserved quantity q can be any component (or linear combination of components) of T^μν and N^μ, and S_q is the corresponding source term. The
explicit equations (<ref>) implemented in  are given in Appendix <ref>. This form allows all quantities to be evolved with the same numerical transport algorithm.

 is designed for flexibility so that different physical limits can be easily studied. Only the propagation of T^ττ, T^τ x, T^τ y, and T^τη 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 ^μ requires the evolution of N^τ to be turned on. 

To solve equation (<ref>) we use the Kurganov-Tadmor (KT) algorithm <cit.>, with a second-order explicit Runge-Kutta (RK) ordinary differential equation solver <cit.> for the time integration step. This scheme is widely used in relativistic hydrodynamic simulations (see, e.g., <cit.>). 




 §.§ 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:


    ∂_τA^n_i,j,k=A^n_i,j,k-A^n-1_i,j,k/Δτ .


Here A is the quantity to be differentiated, i,j,k are  integer labels for the x, y, and η_s coordinates of the grid point, and Δτ is the temporal grid size (numerical resolution in the τ 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^μ is non-zero, for example in the case of the Gubser flow test described in Sec. <ref>.

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


    ∂_xA^n_i,j,k = A^n_i+1,j,k-A^n_i-1,j,k/2Δ x,


where Δ 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 <cit.>.

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:


    ∂_xA^n_i,j,k = (
    θ_fA^n_i,j,k - A^n_i-1,j,k/Δ x, A^n_i+1,j,k - A^n_i-1,j,k/2Δ x, θ_fA^n_i+1,j,k - A^n_i,j,k/Δ x) ,


where the multivariate minmod function is defined as


    (x,y,z)≡(x,(y,z)) ,


with (x,y)≡[(x)+(y)]·(|x|,|y|)/2 and (x)≡ |x|/x. In other words, (x,y,z) always gives the value which is the closest to 0 among (x,y,z).  The parameter θ_f∈[1,2]; θ_f=1 (θ_f=2) corresponds to the most (least) dissipative limiter. In  Eq. (<ref>) is used only for the derivatives of u^μ and , and only when selected by the user as an option.




 §.§ 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. (<ref>) require knowledge of fluid velocity u^μ 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. <cit.>. We here describe the most general form of the root finding problem: assuming that T^τμ, N^τ, π^τμ, Π, and n^τ are all known from the latest temporal update step, we want to compute , ,  and u^μ. As will be demonstrated in Sec. <ref>, the following algorithm <cit.> 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^μ and J^τ to the energy-momentum current T^τμ and baryon density N^τ in the computational frame:


    M^τ   = T^ττ-π^ττ = ( +) ( u^τ )^2 - ,
    
      M^i   = T^τ i-π^τ i = ( +) u^τ u^i   (i = x, y, η_s) ,
    
      J^τ   = N^τ-n^τ = u^τ .


Note that for a viscous fluid =  + Π 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 <cit.> where  is found first, using a 1-dimensional zero search. Here we first find the magnitude of the flow velocity, v, by solving iteratively <cit.>


    v ≡M/M^τ + = 
        M/M^τ +((v),(v)) + Π ,


where M ≡√((M^x)^2+(M^y)^2+τ^2(M^η)^2) and (v), (v) are obtained from the known quantities M^τ, M, and J^τ as


    (v)    = M^τ - vM ,
    
    (v)    = J^τ√(1-v^2) .


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


    u^τ   =1/√(1-v^2) ,
    
    
      u^i    = u^τM^i/M^τ +  ,


for i = x,y,η_s. Note that the algorithm makes active use of all the numerically known components listed above. 

Rewriting Eq. (<ref>) in the form


    f(v) ≡ v-M/M^τ+((v),(v))+Π
        =0 ,


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


    v_i+1 = v_i - f(v_i)/f'(v_i) ,


where


    f'(v) ≡∂ f(v)/∂ v = 1+M/(M^τ+)^2d /d v ,


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


    d /d v=-[M∂/∂+J^τv/√(1-v^2)∂/∂] ,


where (for equations of state like EOS4) the derivatives ∂/∂ and ∂/∂ must be interpolated from the values stored in the EoS table to the pair (,) 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 <cit.> and use Eq. (<ref>) to solve for v only if in the previous time step, at the spatial grid point in question, v ≤ 0.563624 or, equivalently (see (<ref>)), u^τ≤ 1.21061. Otherwise we instead solve for u^τ (which has no upper limit), by employing the Newton-Raphson algorithm to find the zero of


    f(u^τ) ≡ u^τ - √(M^τ+((u^τ),(u^τ))+Π/+((u^τ),(u^τ))+Π) .


In this case we need to evaluate in each iteration 


    f'(u^τ) = 1 - 1/2[-M^τ/(+)^3/2(M^τ+)^1/2]d/du^τ ,


where


    d/du^τ = -1/(u^τ)^2[M/vu^τ∂/∂ + J^τ∂/∂] ,


with v = √(1-1/(u^τ)^2). In Sec. <ref> it will be shown that the switch between the two schemes, as implemented in , works seamlessly, moving smoothly from v = 0.563624 to u^τ = 1.21061 across the switching point. 

When interpolating the EoS table to obtain the derivatives needed on the r.h.s. of Eqs. (<ref>,<ref>) one can encounter numerical errors in regions where derivatives of the EoS change discontinuously, e.g. in the recently developed BEST EoS <cit.> which adds (using a certain prescription) a critical point and first-order phase transition to the LQCD-HRG interpolated EOS4. For such situations,  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 ((v_i), (v_i)) from Eqs. (<ref>,<ref>) and the EoS (,), compute an updated value v_i+1 of the velocity from


    v_i+1 = M/M^τ +((v_i),(v_i)) + Π ,


and iterate these steps until convergence is reached. For v ≥ 0.563624 or u^τ≥ 1.21061, one instead updates u^τ using the equation 


    u^τ_i+1 = 
      √(M^τ+((u^τ_i),(u^τ_i))
            +Π/ +((u^τ_i),(u^τ_i))+Π)


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. <ref> 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 π^xx, π^xy, π^xη, and ^η. In principle, these could be computed from the other components of the shear stress and baryon diffusion current by using the tracelessness of π^μν and the orthogonality of π^μν and ^μ to u^μ. 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.




 §.§ 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 <cit.>, this may not happen quickly enough to avoid numerical instability of the evolution algorithm.[
    Note 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 <cit.>.]

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) <cit.>. 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.[
    For example, CLVisc <cit.> requires (|π^μν|) < T_0^ττ; when this is violated for some cell in the dilute region, π^μν is set to 0 locally. vHLLE <cit.> requires (|π^μν|/|T_0^μν|) < C and |Π|/ < C, with C being a constant of order but smaller than 1; if one of these conditions is violated, Π and/or π^μν are rescaled by a factor (which is common for all components of π^μν) to satisfy this requirement.]

In , we follow the lead of iEBE-VISHNU <cit.>, GPU-VH <cit.> and MUSIC <cit.> 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) π^μν with T_0^μν= u^μ u^ν - Δ^μν, Π with √(^2+3^2), and ^μ with u^μ. 

For the shear stress tensor, iEBE-VISHNU <cit.> and GPU-VH <cit.> require 


    √(π^μνπ_μν)≤ρ_max√(T_0^μνT_0μν) = ρ_max√(^2+3^2) ,


with ρ_max≤ 1. In addition, the tracelessness of π^μν and its orthogonality to u_ν are required,


    π^μ_μ≤ξ_0√(π^μνπ_μν)  and  π^μνu_ν≤ξ_0√(π^μνπ_μν) ,


where √(π^μνπ_μν) sets the scale and ξ_0≪ 1 is a small number <cit.>. At grid points where these trigger conditions are violated, Refs. <cit.> regulate the shear stress tensor π^μν by (see the left plot of Fig. <ref>)


    π^μν→tanhρ_π/ρ_ππ^μν ,


where


    ρ_π≡[
      √(π^μνπ_μν)/ρ_max√(^2+3^2) ,  g_μνπ^μν/ξ_0ρ_max√(π^μνπ_μν) ,  π^λμu_μ/ξ_0ρ_max√(π^μνπ_μν) ∀ λ].


For the bulk viscous pressure, which can make the root finding process fail when negative and too large, Refs. <cit.> regulate Π during the root finding process to ensure existence of at least one non-negative solution for , v or u^τ. In  we demand instead that 


    √(3Π^2)≤ρ_max√(^2+3^2) ,


and where this trigger condition is violated we regulate Π by


    Π→tanhρ_Π/ρ_ΠΠ ,


with


    ρ_Π≡1/ρ_max√(3Π^2/^2+3^2) .



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


    -^μ_μ≪^2   and  ^μ u_μ = 0 .


In the code we replace these conditions by

    √(-^μ_μ)≤ρ_max√(^2)  and  ^μ u_μ≤ξ_0√(-^μ_μ) .

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


    ^μ→tanhρ_/ρ_^μ ,


with


    ρ_≡[
      √(-^μ_μ)/ρ_max√(^2) ,  ^μ u_μ/ξ_0ρ_max√(-^μ_μ)] .







Equations (<ref>,<ref>,<ref>) define the first of our two regulation schemes. In contrast to Refs. <cit.> where Π 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 Π during the root finding process. As defaults for the regulation parameters we use the same values ξ_0 = 0.1 and ρ_max = 1 as proposed in <cit.>. For documentation of the sensitivity studies leading to these default values we refer the reader to Ref. <cit.>.

A second, different regulation scheme is based on the one implemented in the latest version of MUSIC <cit.>. (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


    π^μν→r^max_π/r_ππ^μν ,  Π→r^max_Π/r_ΠΠ ,  ^μ→r^max_/r_^μ ,


where regulation of quantity i is triggered whenever r_i exceeds the corresponding maximally allowed value r_i^max (i=π, Π, ),[
    If the regulation for the shear stress or baryon diffusion is triggered, all components of π^μν or n^μ are regulated by a common regulation factor.]

with r_i defined by 


    r_π = 1/f_s√(π^μνπ_μν/^2+3^2) ,   
        r_Π = 1/f_s√(3Π^2/^2+3^2) ,   
        r_ = 1/f_s√(-^μ_μ/^2) .


By comparing Eq. (<ref>) to Eqs. (<ref>,<ref>,<ref>) and Eq. (<ref>) to Eqs. (<ref>,<ref>,<ref>) above, we see that the factors r_i^max/r_i play the same role as tanhρ_i/ρ_i in the first regulation scheme, causing stronger regulation for larger r_i, with r_i playing the role of the quantity ρ_i whereas f_s, defined by


    f_s = χ_0
        [1/exp(-(-_0)/ξ_0)+1
             -1/exp(_0/ξ_0)+1],


playing a similar role as ρ_max: as f_s or ρ_max grows larger, the regulation gets weaker. f_s is designed to approach χ_0 when ≫_0 and 0 when ≪_0.

The right plot of Fig. <ref> shows that for < _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 χ_0, the regulation will hardly ever be triggered in the dense region >_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 _0=0.1 GeV/fm^3 <cit.> 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. <cit.> used χ_0=10, and chose r^max_ = 1 for regulating the baryon diffusion current; we adopt the same value for χ_0 and identical maximum r values for all dissipative flows: r^max_π = r^max_Π = r^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. <cit.>).




§ 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 <cit.> was carefully validated using similar tests, but the entire baryon evolution sector in  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 <cit.> also includes baryon evolution) will become possible when that version of MUSIC becomes public.

Building on validation protocols described in Refs. <cit.>, we here discuss tests in which we compare, for identical initial conditions, the numerical solutions from  with (semi-)analytic solutions using Mathematica <cit.> for the Riemann problem for the Euler equations <cit.>, for Bjorken flow <cit.>, and for Gubser flow <cit.> 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  with the independently developed numerical algorithm described in Ref. <cit.> 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. <ref> 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 ħ 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: Δ x = Δ y = 0.05 fm and Δη_s = 0.02; Δτ 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.




 §.§ 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 <cit.>. Here a special case of the relativistic Riemann problem, known as Sod's shock tube problem <cit.>, 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:


    (,,u^x) =
        
             (_l, _l, u^x_l = 0) ,     x < 0 ,
    
             (_r, _r, u^x_r = 0) ,     x > 0 .


In the longitudinal (z) direction the fluid is assumed to expand with a boost-invariant velocity profile u^z/u^τ = z/t (i.e. u^η = 0), and the system is initialized along a surface of constant longitudinal proper time τ_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 <cit.> (see Fig. <ref>(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. <ref>(a,c)) while the density has a jump (Fig. <ref>(b)) <cit.>. 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., <cit.>). Since for EOS1 ∂/∂ (which is needed in Eqs. (<ref>,<ref>) 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  and T are obtained from the EoS. For ideal fluid dynamics the equations of motion become


    D =    - θ , 
    
        D =    -(+)θ ,
    
        (+)Du^μ =    - Δ^μν∂_ν .


They are solved with initial conditions (<ref>) with the default parameters listed in Fig. <ref>. For details about the analytical solution we refer the reader to Ref. <cit.> (see also <cit.>). We point out that in the analytical solution from Ref. <cit.> the evolution starts at time zero whereas in the code the hydrodynamic evolution is initialized at τ_0 = 0.5 fm. Since the solution is self-similar and depends only on the variable x/(τ-τ_0), we therefore compare in Fig. <ref> the numerical results at τ=8.5 fm/c to the analytical solution at τ-τ_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. (<ref>), and not Eq. (<ref>) which would yield zero initial gradients and result in no evolution at all.

Figure <ref> 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. 




 §.§ Bjorken flow




In this subsection we test  in Milne coordinates for a transversally homogeneous dissipative fluid undergoing longitudinally boost-invariant Bjorken expansion <cit.>. Boost-invariant systems are characterized by space-time rapidity independent macroscopic observables and a flow profile that looks static (i.e. u^μ = (1,0,0,0)) in Milne coordinates <cit.>. 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 √(s)≳ 100 GeV/nucleon (see, e.g., <cit.>). 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:


    =   -++Π -π/τ ,
    
    τ _ΠΠ̇+Π   =   -ζ/τ-δ _ΠΠ
    Π/τ+λ _Πππ/τ ,
    
    τ _ππ̇+π   =   4/3η/τ -( 1/3τ
    _ππ+δ _ππ) π/τ+2/3λ _πΠ
     Π/τ ,
    
    
    
    
       =   -/τ ,
    
    τ_^η+^η   =   -(τ_+δ_+2/3λ_)^η/τ ,


where π≡-τ^2π^ηη has been introduced. 






In the  simulation the numbers of cells in (x, y, η_s) directions are all set to 1.[
    Of course, one can also set the number of cells in (x, y) larger than 1 but to ensure transverse spatial homogeneity quantities like  and  should be the same in all cells.]

With boost-invariant and transversally homogeneous initial conditions, the numerical results from the (3+1)D  code can be tested against a solution of the coupled ODEs (<ref>-<ref>) from a separate Mathematica code. As for the Riemann problem, the baryon evolution decouples from the rest of the system if an EoS = () 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. (<ref>,<ref>), without affecting the Bjorken flow profile.  

By using EOS2 (the Wuppertal-Budapest EoS at μ=0) which features a non-zero interaction measure, we can also test the evolution of the bulk viscous pressure Π which is propagated by Eq. (<ref>) as an additional dissipative degree of freedom <cit.>. Noting that Eq. (<ref>)
involves a rather complex parametrization of the transport coefficients (<ref>)-(<ref>), we emphasize that this test is indeed non-trivial.

As Fig. <ref> shows, the agreement of  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.[
    As for the Riemann problem, this test again checks this property only for the case ∂/∂ = 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”.





 §.§ Gubser flow





Gubser flow <cit.> 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 τ_0). 




  §.§.§ 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⊗ℝ <cit.>. 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, ρ∈ℝ.

To introduce the Gubser coordinates we first rescale the invariant distance measure of Minkowski space in Milne coordinates with a Weyl transformation[
    To 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^μν = diag(-1,1,1,1) in Cartesian coordinates.]  


    ds^2→ dŝ^2≡ ds^2/τ^2 = (-dτ^2+dr^2+r^2dϕ^2)/τ^2+dη_s^2 .


Next we perform the coordinate transformation x^μ = (τ, r, ϕ, η_s) →x̂^μ = (ρ,θ,ϕ, η_s),[
    Here r≡√(x^2+y^2) and ϕ≡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 τ and labeled with a hat.] 

by introducing <cit.>


    ρ(τ,r)    ≡ - sinh^-1(1-q^2τ^2+q^2r^2/2qτ) ,
    θ(τ,r)    ≡tanh^-1(2qr/1+q^2τ^2-q^2r^2) ,

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, τ and r such that qr and qτ remain unchanged). In these coordinates the Weyl-rescaled invariant distance measure becomes


    dŝ^2 = -dρ^2+cosh^2ρ (dθ^2+sin^2θ dϕ^2)+dη_s^2 ,


with the metric 


    ĝ_μν = diag(-1, cosh^2ρ, cosh^2ρsin^2θ,  1) .


A system that appears static in the coordinates x̂^μ = (ρ, θ, ϕ, η_s), i.e. has flow velocity û^μ=(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 <cit.> and the definitions (<ref>,<ref>), for example


    u_μ(τ, r)    =τ∂x̂^ν/∂ x^μû_ν(ρ(τ, r)) ,
    
    π_μν(τ, r)    =1/τ^2∂x̂^α/∂ x^μ∂x̂^β/∂ x^νπ̂_αβ(ρ(τ, r)) ,
    
    
    T(τ, r)    = T̂(ρ(τ, r))/τ ,
    
    (τ, r)    = (ρ(τ, r))/τ^4  .


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


    u^τ(τ, r)    = coshκ(τ,r) , 
    
    
    u^x(τ, r)    = x/rsinhκ(τ,r) ,
    
    
    u^y(τ, r)    = y/rsinhκ(τ,r) ,
    
    
    u^ϕ(τ, r)    = u^η(τ, r) = 0 ,


where κ(τ,r) is the transverse flow rapidity, corresponding to the transverse flow velocity 


    v_⊥(τ,r)=tanhκ(τ,r) ≡2q^2τ r/1+q^2τ^2+q^2r^2 .


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 ρ when expressed in Gubser coordinates). Different initial conditions for the hydrodynamic quantities and different transport coefficients yield different ρ dependencies for their evolution, translating into different characteristics r(τ) and different (τ, r(τ)) profiles when expressed in Minkowski space coordinates. 

To ensure invariance of the hydrodynamic equations under the Weyl transformation (<ref>) 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.
 



  §.§.§ Gubser flow with baryon diffusion




In this work we extend the existing semi-analytical solutions for conformal Israel-Stewart hydrodynamics with Gubser flow <cit.> to systems with a non-zero baryon diffusion current. One could argue that the longitudinal reflection symmetry under η_s→-η_s in the Gubser symmetry indicates that the baryon diffusion current should be zero <cit.>. 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 û^μ=(1,0,0,0) under Gubser symmetry. 

In this subsection we derive these baryon equations of motion in de Sitter space. We rewrite Eqs. (<ref>, <ref>) and (<ref>, <ref>) with the mostly-plus metric tensor and apply the Gubser flow profile û^μ=(1,0,0,0) and the de Sitter metric (<ref>) to obtain[
    Similar to Bjorken flow, the shear stress for Gubser flow has only one independent component for which we take π̂^ηη. The other non-vanishing components π̂^θθ and π̂^ϕϕ are related to π̂^ηη by tracelessness (which gives π̂^η_η = -π̂^θ_θ - π̂^ϕ_ϕ) and azimuthal symmetry (which implies π̂^η_η = -2π̂^θ_θ = -2π̂^ϕ_ϕ for evolution with azimuthally symmetric initial conditions).
    ]


    ∂_ρ+ 2tanhρ    =    
       2tanhρ (1/2π̂^ηη - P̂_0),
    
    τ _π∂_ρπ̂^ηη+π̂^ηη   =    2tanhρ(2/3η̂-δ_πππ̂^ηη+1/6τ_πππ̂^ηη),  
    
    ∂_ρ+ 2tanhρ    =   0 ,
    
    τ _∂_ρ^η +^η   =    -2tanhρ(δ_ - 1/3λ_) ^η .

(Note that θ̂=2tanhρ is the scalar expansion rate for Gubser flow.) For the transport coefficients we use the same parametrization as described before in Fig. <ref> (see also the caption of Fig. <ref>). The transformation rules for the baryon density and diffusion current are


    (τ, r)    =   (ρ(τ, r))/τ^3  ,
    
    ^η(τ, r)    =   ^η(ρ(τ, r))/τ^4 .







Given initial conditions (ρ_0), π̂^ηη(ρ_0), (ρ_0), and ^η(ρ_0) at, say, ρ_0 = 0, obtaining (ρ_0) and T̂(ρ_0) with the help of EOS1, Eqs. (<ref>-<ref>) can be solved with Mathematica in Gubser coordinates and then translated into Milne coordinates (τ,x,y,η_s=0) using the transformation rules (<ref>-<ref>) and (<ref>-<ref>). This semi-analytic solution can then be used to obtain initial conditions for  on an initial proper time hypersurface τ_0 which are then further evolved with the (3+1)-dimensional  code. Specifically,  requires initial data at τ_0 for 
    
    , , T, u^x, u^y, π^xx, π^yy, π^xy, π^ττ, π^τ x, π^τ y, π^ηη, , n^η

on the computational (x, y) grid (due to longitudinal boost-invariance these are only required at η_s=0). All remaining hydrodynamic components are either zero or can be obtained from the above by symmetry. For example, all shear stress components π^μν(τ_0, x, y, η_s=0) for  can be obtained from the semi-analytic solution π̂^ηη(ρ) by using tracelessness and azimuthal symmetry,


    π̂_θθ(ρ)   =-1/2cosh^2ρ π̂_ηη(ρ) ,   
        π̂_ϕϕ(ρ)=-1/2cosh^2ρsin^2θ π̂_ηη(ρ)


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


    π_μν(ρ(τ, r(x, y))) =1/τ^2∂x̂^α/∂ x^μ∂x̂^β/∂ x^νπ̂_αβ(ρ) .

This gives, for example, at (τ=τ_0,η_s=0)

    π^ττ(ρ(τ_0, r(x, y))) = -q^2sin^2θ/2τ_0^2 π̂_ηη(ρ) ,

where the value of ρ depends on the transverse grid point (x,y).[
    Obviously, whenever the transport coefficients are changed in , 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  code in Milne coordinates.]














In Fig. <ref> we compare the  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  evolution equations when compared with the Bjorken flow test. We can also use it to study the performance of the root-finding algorithm in . 

As discussed in Sec. <ref>, two different methods are used in complementary ranges of the flow velocity separated by the critical value v = 0.563624 or, equivalently, u^τ = 1.21061. Fig. <ref>(e) shows that the root-finding algorithm works equally well on both sides of the critical value of u^τ. We also checked the precision and relative speed of convergence of the Newton-Raphson and modified iteration schemes described in Sec. <ref>. 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^τ_i+1-u^τ_i|/u^τ_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.[
    During 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.] 







 §.§ 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  code in describing the expansion along the rapidity direction. To remedy this we have compared  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 <cit.>. To be able to study baryon number transport, EOS4 is used in both codes. A similar comparison was also made in Ref. <cit.> to test the performance of MUSIC. Rather than directly using Monnai's code <cit.> we compare our  results with those reported in the comparison <cit.> with MUSIC. In this sense, the following test is also a code comparison with MUSIC. Qualitatively consistent results from an earlier  study were already reported in <cit.>.






We here focus on the baryon number evolution in the longitudinal direction, with non-vanishing longitudinal gradients of μ/T, by setting bulk and shear stresses to zero.[
    In 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  handles the evolution equation (<ref>) for baryon number diffusion,


    d ^μ = κ_n/τ_∇^μ(μ/T) - ^μ/τ_ - n^ν u^μ D u_ν 
        - u^αΓ^μ_αβ n^β ,

together with the net baryon conservation law (<ref>). For the baryon transport coefficients we choose κ_n = 0.2 /μ and τ_n = 0.2/T. 

In Fig. <ref> we show a comparison of the distributions in space-time rapidity η_s of energy density  and net baryon density  at four different times. Broken (continuous) lines show the
results from  (MUSIC <cit.>); 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  (,) that depends on both energy and net baryon density. 




 §.§ Tests summary



We briefly summarize which parts of the  code were tested with the test protocols described in this section. As described in Sec. <ref>, 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. <ref>) shows the ability of the RK-KT algorithm in  to capture shocks and contact discontinuities. The (0+1)D Bjorken expansion with EOS2 (Sec. <ref>) 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. <ref>) 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. <ref>) 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  passes all these tests without struggle.




§ 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 η_s=0.[
    For a similar earlier study with smooth, ensemble-averaged initial conditions see <cit.>.]
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. <cit.>), we here use smooth longitudinal initial conditions and refrain from showing the (mostly uninteresting) evolution along the beam direction.


Following Refs. <cit.> 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 <cit.> into the longitudinal direction with the following prescription:


    (τ_0,x,y,η_s)    =   _0/τ_0 [T_A(x,y)_A(η_s)+T_B(x,y)_B(η_s)] ,
    
    (τ_0,x,y,η_s)    =   1/τ_0 [T_A(x,y)_A(η_s)+T_B(x,y)_B(η_s)] ,


where T_A/B(x,y) are the transverse profiles of the right- and left-moving nuclei from the MC-Glauber model and _A/B(η_s), _A/B(η_s) are the corresponding longitudinal profiles for the energy and net baryon density, respectively <cit.>. _0 is a normalization factor which can be tuned to reproduce the final multiplicity while (τ_0,x,y,η_s) is normalized to the total number of participant baryons <cit.>.





In Fig. <ref>, we show the transverse distributions at η_s=0 of the energy density (upper panels) and net baryon density (lower panels), at the hydrodynamic starting time τ_0 =0.5 fm/c (left column) and later at τ=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. <cit.> came to a similar conclusion for the evolution in the longitudinal direction.




§ CONCLUSIONS AND OUTLOOK



In this work we have described the structure and performance of a new code called  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.  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  algorithm to the simultaneous propagation of multiple conserved charge currents <cit.> should be a straightforward task for the near future. A dynamical initialization module for  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. <cit.>, is under construction and will soon be added to the  package. The code can also be plugged in the JETSCAPE framework <cit.> as a hydrodynamic evolution module.

After briefly describing the physics to be addressed with  simulations, we presented in detail the equations  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 <cit.> and by the Open Science Grid <cit.>, which is supported by the National Science Foundation, award 1148698, and the U.S. Department of Energy's Office of Science.







§ CONSERVATIVE FORM OF THE EVOLUTION EQUATIONS



In this Appendix, we recast Eqs. (<ref>-<ref>) and  (<ref>-<ref>) in the conservative form (<ref>). The procedure follows Refs. <cit.>, adding here the equations for baryon evolution. The equations reproduced here are written in a form that facilitates direct comparison with the  code.

We start with the conservation laws in Eqs. (<ref>-<ref>). With the scaled flow velocities v^i=u^i/u^τ (i=x,y,η_s), we can write down the following constituent relations for the components of T^μν and N^μ:


    T^ττ   =(+)u^τu^τ-+π^ττ ,
    
    
    T^τ i   =(+)u^τu^i+π^τ i =v^iT^ττ+ v^i-v^iπ^ττ+π^τ i ,
    
    
    T^ij   =(+)u^iu^j- g^ij+π^ij =v^iT^τ j- g^ij-v^iπ^τ i+π^ij ,
    
    
    
    
    N^τ   =  u^τ + n^τ , 
    
    N^i    = u^i + n^i = v^i N^τ - v^i n^τ + n^i .

Inserting these into the conservation laws one obtains


    ∂ _τT^ττ+∂ _x(v^xT^ττ)+∂_y(v^yT^ττ)+∂ _η(v^ηT^ττ)    = I^τ_2 + I^τ_x +I^τ_y +I^τ_η , 
    
    
    ∂ _τT^τ x+∂ _x(v^xT^τ x)+∂
    _y(v^yT^τ x)+∂ _η(v^ηT^τ x)    = I^x_2 + I^x_x +I^x_y +I^x_η ,
    
    
    ∂ _τT^τ y+∂ _x(v^xT^τ y)+∂
    _y(v^yT^τ y)+∂ _η(v^ηT^τ y)    = I^y_2 + I^y_x +I^y_y +I^y_η , 
    
    
    ∂ _τT^τη+∂ _x(v^xT^τη)+∂ _y(v^yT^τη)+∂ _η(v^ηT^τη)    = I^η_2 + I^η_x +I^η_y +I^η_η , 
    
    ∂ _τN^τ+∂ _x(v^xN^τ)+∂ _y(v^yN^τ)+∂ _η(v^ηN^τ)   = J^τ_2 + J^τ_x +J^τ_y +J^τ_η ,

with the following source terms for T^ττ:

    I^τ_2=   -1/τ( T^ττ+τ ^2T^ηη)-(+Π-π ^ττ)∂ _iv^i-v^i∂_i  , 
    
    I^τ_x=    -v^x∂ _x(Π-π ^ττ)-∂ _xπ ^τ x , 
    
    I^τ_y=    -v^y∂ _y(Π-π ^ττ)-∂ _yπ ^τ y , 
    
    I^τ_η=    -v^η∂ _η(Π-π ^ττ)-∂ _ηπ ^τη ;


for T^τ x:

    I^x_2=   -1/τT^τ x-∂_x +π ^τ x∂ _iv^i , 
    
    I^x_x=    -∂ _x( Π+π ^xx)+v^x∂ _xπ ^τ x ,
    
    I^x_y =      -∂_yπ ^xy+v^y∂_yπ ^τ x , 
    
    I^x_η=     -∂_ηπ ^xη+v^η∂_ηπ ^τ x ;


for T^τ y:

    I^y_2=   -1/τT^τ y-∂_y +π ^τ y∂ _iv^i ,
    
    I^y_x=   -∂ _xπ ^xy+v^x∂ _xπ ^τ y ,
    
    I^y_y=   -∂_y(Π+π ^yy)+v^y∂_yπ ^τ y , 
    
    I^y_η=   -∂_ηπ ^yη+v^η∂_ηπ ^τ y ;


for T^τη:

    I^η_2=   -3/τT^τη
              -∂_η/τ^2
              +π^τη∂_iv^i ,
    
    
    I^η_x=   -∂_xπ^xη+v^x∂_xπ^τη ,
    
    
    I^η_y=   -∂_yπ ^yη+v^y∂_yπ^τη ,
    
    
    I^η_η=   -∂_η(Π/τ ^2+π^ηη)
    + v^η∂_ηπ^τη ;


and for N^τ:

    J^τ_2=   -1/τ N^τ+n^τ∂ _iv^i  ,
    
    
    J^τ_x=   -∂ _xn ^x +v^x∂ _xn ^τ ,
    
    
    J^τ_y=    -∂ _yn^y +v^y∂ _yn^τ ,
    
    
    J^τ_η=    -∂_ηn^η 
               +v^η∂_ηn^τ .


Considering d ≡ u^μ∂_μ, the relaxation equations (<ref>-<ref>) for the dissipative flows can be written as

    ∂_τΠ+∂_x(v^xΠ)+∂_y(v^yΠ)+∂_η(v^ηΠ)   =   S^Π_2  ,
    
    ∂_τ^μ+∂_x(v^x^μ)+∂_y(v^y^μ)+∂_η(v^η^μ)    =   S^_2
     ,
    ∂_τπ^μν +∂_x(v^xπ^μν)+∂_y(v^yπ^μν)+∂_η(v^ηπ^μν)    =   S^π_2  ,

where we used ∂_iv^i≡∂_xv^x + ∂_yv^y + ∂_ηv^η. The source terms in Eqs. (<ref>-<ref>) are given by 


    S^Π_2   =   1/u^τ(-ζ/τ_Πθ-Π/τ_Π-I_Π)+Π∂_i v^i
         ,
    
        S^_2   =   1/u^τ(κ_n/τ_∇^μ(μ/T)-^μ/τ _-I_^μ-G_^μ)+^μ∂_i v^i 
         ,
    
        S^π_2   =   1/u^τ(2η/τ_πσ^μν-π^μν/τ_π-I^μν_π-G^μν_π)+π^μν∂_i v^i .

Here the I-terms and G-terms are defined in Eqs. (<ref>)-(<ref>).




§ EXPLICIT FORM



In this Appendix we provide explicit forms of some expression occurring in the baryon evolution equations. In Milne coordinates the terms G^μ_n = u^αΓ^μ_αβ n^β in Eq. (<ref>) evaluate to 


    G^τ_n    = τ u^η n^η ,
    
    G^x_n    = 0 ,
    
    G^y_n    = 0 ,
    
    G^τ_n    = (u^τ n^η + u^η n^τ)/τ .

To calculate the LRF gradient of μ/T in Eq. (<ref>) numerically we use


    ∇^μ(μ/T) ≡Δ^μν d_ν(μ/T) = (g^μν-u^μ u^ν)∂_ν(μ/T)

and work out the partial derivatives ∂_μ (μ/T) in the computational frame numerically from the EoS tables. Finally, the last two terms in Eq. (<ref>) can be expressed as


    I_3^μ    = n_νω ^νμ = n^τω ^τμ - n^xω ^x μ - n^yω ^y μ - τ^2 n^ηω ^ημ ,
    
    I_4^μ    = n_νσ ^νμ= n^τσ^τμ - n^xσ^x μ - n^yσ^y μ  - τ^2 n^ησ^ημ .








elsarticle-num

