Baryon transport and the QCD critical point

Fireballs created in relativistic heavy-ion collisions at different beam energies have been argued to follow different trajectories in the QCD phase diagram in which the QCD critical point serves as a landmark. Using a (1+1)-dimensional model setting with transverse homogeneity, we study the complexities introduced by the fact that the evolution history of each fireball cannot be characterized by a single trajectory but rather covers an entire swath of the phase diagram, with the finally emitted hadron spectra integrating over contributions from many different trajectories. Studying the phase diagram trajectories of fluid cells at different space-time rapidities, we explore how baryon diffusion shuffles them around, and how they are affected by critical dynamics near the QCD critical point. We find a striking insensitivity of baryon diffusion to critical effects. Its origins are analyzed and possible implications discussed.


I. INTRODUCTION
One of the primary goals of nuclear physics [1] is studying the phase diagram of Quantum Chromodynamics (QCD), which is generally mapped to a plane expanded with temperature T and baryon chemical potential µ axes [2]. First principles calculations from Lattice QCD show that, at zero µ, the phase transition from a deconfined quark-gluon plasma (QGP) phase to a confined hadron resonance gas (HRG) phase from high to low temperature is a rapid but smooth crossover [3][4][5][6]. At large µ, calculations of the phase transition using Lattice QCD are not yet available, since there the standard techniques suffer from the "sign problem" [7,8]. Nevertheless, theoretical models indicate that at large chemical potential, the phase transition is first order [2], and this implies that a critical point exists at non-zero chemical potential [9,10], at the end of the first order phase transition line. Confirming the existence and finding the location of the hypothetical QCD critical point have attracted tremendous amount of attention over the last two decades [11,12].
Heavy-ion collisions are the main method to tackle these unsolved problems [9][10][11][12][13][14][15][16][17][18][19][20][21][22]. Such collisions have been carried out at different experimental facilities, such as the Large Hadron Collider (LHC) at CERN and the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Laboratory, at various beam energies, and large sets of data have been accumulated. One of the most promising signatures of the QCD critical point is a nonmonotonic beam energy dependence of higher-order cumulants of the fluctuations in the net proton production yields [9,10,16,17,19,20]. This is based on the idea that these observables are more sensitive to the correlation length of fluctuations of the chiral order parameter which, in the thermodynamic limit, diverges at the critical point [23]. Fireballs created in heavy-ion collisions at different beam energies should freeze out with correlation lengths that depend non-monotonically on the collision energy, and this should be reflected in the net baryon cumulants. Strongly motivated by this, a Beam Energy Scan (BES) program has been carried out at RHIC during the last decade. During a first campaign that ended in 2011 (BES-I), Au-Au collisions were studied at collision energies √ s NN from 200 GeV down to 7.7 GeV (BES-I) [12,24,25]. A second campaign, BES-II, with significantly increased beam luminosity is expected to be completed this year, after having explored collision energies down to √ s NN = 3.0 GeV in fixed-target mode. Additional experiments at even lower beam energies, probing the phase diagram in regions with even higher baryon chemical potential, are planned at the newly constructed FAIR and NICA facilities. Unfortunately, the dynamical nature of the fireballs created in heavy-ion collisions renders attempts to confirm the above static equilibrium considerations experimentally anything but straightforward. Within their short lifetimes of several dozen yoctoseconds the fireballs' energy density decreases rapidly by collective expansion, from initially hundreds of GeV/fm 3 to well below 1 GeV/fm 3 at final freeze-out (see e.g. Refs. [11,26]). The rapid dynamical evolution of the thermodynamic environment keeps the latter permanently out of thermal equilibrium such that critical fluctuations never reach their thermodynamic equilibrium distributions. In addition, in those parts of the fireball which pass through the quark-hadron phase transition close to the QCD critical point the dynamics of critical fluctuations is affected by "critical slowing down" [27]. This is both a curse and a blessing: If critical fluctuations would relax quickly to thermal equilibrium, all memory of critical dynamics might have been erased from the hadronic freeze-out distributions by the time the hadron yields and momenta decouple. If, on the other hand, the dynamical evolution of fluctuations is slowed in the vicinity of the critical point, some signals of critical dynamics may survive until freeze-out but they will then most definitely not feature their equilibrium characteristics near the critical point [27].
Thus, to confirm or exclude the critical point via systematic model-data comparison, reliable dynamical simulations of off-equilibrium critical fluctuations and the associated final particle cumulants, on top of a wellconstrained comprehensive dynamical description of the bulk medium at various beam energies, are indispensable [28][29][30][31][32][33][34][35][36][37][38][39]. Recently, the Hydro+/++ framework [34,40] was developed for incorporating off-equilibrium fluctuations and critical slowing-down into hydrodynamic simulations, and some practical progress within simplified settings has since been made using this framework [21,22]. On the other hand, while a fully developed and calibrated (2+1)-dimensional multi-stage description of heavy-ion collisions (including initial conditions + pre-hydrodynamic dynamics + viscous hydrodynamics + hadronic afterburner) exists (see, e.g., the most recent versions described in [41][42][43][44]) and has met great phenomenological success at top RHIC and LHC energies [11,12,26], such a comprehensive and fully validated framework is still missing for collisions at the lower BES energies.
Compared to high-energy collisions at LHC and top RHIC energies, collisions at BES energies introduce a number of additional complications [45,46]. These include (i) a much more complex, intrinsically (3+1)dimensional and temporally extended nuclear interpenetration stage and its associated dynamical deposition of energy and baryon number [47][48][49][50][51], (ii) the need to account for and properly propagate conserved charge currents for baryon number and strangeness [51][52][53][54][55][56][57], (iii) a consistent treatment of singularities in the thermodynamic properties associated with the critical point [23,58], (iv) the aforementioned off-equilibrium nature of critical dynamics [21,22,34,38,40,59], and (v) the dynamical effects of nucleation and spinodal decomposition in the metastable and unstable regions associated with the first-order phase transition [60][61][62][63]. The situation is made even more complicated by the back-reaction of the non-equilibrium critical fluctuation dynamics on the bulk evolution of the medium. This back-reaction causes a potential dilemma: On the one hand, locating the critical point requires reliable calculations at various beam energies of critical fluctuations on top of a wellconstrained bulk evolution of the fireball medium; on the other hand, the back-reaction of the off-equilibrium critical fluctuations from a critical point whose location is yet to be determined onto the medium evolution might interfere with the calibration of the latter and turn it into an impossibly complex iterative procedure whose convergence cannot be guaranteed.
Guidance on how (and perhaps even whether) to incorporate critical effects when constraining the bulk dynamics is direly needed, not least since dynamical simulations for low beam energies are computationally very expensive. Some effects on the bulk medium evolution arising from singularities in the thermodynamic properties of the QCD matter [23,58], by adding a critical point to the QCD Equation of State (EoS) [53,[64][65][66][67] and/or explicitly including critical scaling of its transport coefficients, have been explored, with special attention to the bulk viscous pressure since critical fluctuations of the chiral order parameter, which couples to the baryon mass, can be directly related to a peak of the bulk viscosity near the critical point [68][69][70][71]. The authors of [72] showed that critical effects on the bulk viscous pressure have non-negligible phenomenological consequences for the rapidity distributions of hadronic particle yields, implying that critical effects might indeed play an important role in the calibration of the bulk medium.
In a similar spirit we study here critical effects on the bulk evolution in the baryon sector, by including the critical scaling of the relaxation time for the baryon diffusion current and of the baryon diffusion coefficient, as well as (in a simplified treatment) the critical contribution to the EoS [64,65,67]. We study the phenomenological consequences of baryon diffusion in a system both away from and close to the QCD critical point; the former is essential for modeling heavy-ion collisions at the high end of the BES collision energy range [51][52][53][54][55][56][57]. Including only the baryon diffusion while neglecting other dissipative effects helps us to study its hydrodynamical consequences in isolation. A more comprehensive study including all dissipative effects simultaneously is left for future work.
This paper is organized as follows. We discuss the hydrodynamic formalism with non-zero baryon diffusion current in Sec. II. In Sec. III, we illustrate our setup near the QCD critical point, particularly by discussing the critical behavior of thermodynamic and transport coefficients for the hydrodynamic evolution. After completing the setup of the framework in Sec. IV we discuss results for a fireball created at low beam energies in Sec. V, with a focus on baryon diffusion current effects, first away from (Sec. V A) and then close to (Sec. V C) the critical point, followed by a discussion of general features of the timeevolution of the diffusion current in Sec. V D. We summarize results and draw conclusions in Sec. VI. In the Appendices, we discuss causality near the QCD critical point in App. A, estimate the size of the critical region in App. B and, finally, validate the numerical methods used in this work in App. C.
Throughout this article we use natural units where factors ofh, c, and k B are not explicitly exhibited but implied by dimensional analysis. We also use Milne coordinates, x µ = (τ, x, y, η s ), where τ and η s are the (longitudinal) proper time and space-time rapidity, respectively, and are related to the Cartesian coordinates via t = τ cosh η s , z = τ sinh η s . We employ the mostly-minus convention for the metric tensor, g µν = diag (+1, −1, −1, −1/τ 2 ).

II. VISCOUS HYDRODYNAMICS WITH BARYON DIFFUSION
In this section we discuss the viscous hydrodynamic framework including baryon diffusion. Hydrodynamics is an effective theory for describing long-wavelength degrees of freedom, which, from a macroscopic viewpoint, can be formulated by conservation laws of hydrodynamic variables that are ensemble-averaged at certain coarsegrained scales [26,73]. In heavy-ion collisions, the conserved quantities are energy, momentum, and various charges, including net baryon charge, electric charge and strangeness. In this work, we only study the net baryon charge; the extension to incorporate other conserved charges is conceptually straightforward [56] and will be studied elsewhere. The conservation equations for energy-momentum and net baryon charge are formulated covariantly as where d µ is the covariant derivative in an arbitrary coordinate system (Milne coordinates here), T µν and N µ are the energy-momentum tensor and (net) baryon current, respectively. In a given arbitrary reference frame, they can be decomposed into the ideal and dissipative parts as where T µν id and N µ id are the ideal parts which are well defined in local equilibrium, while Π µν and n µ are the dissipative components describing the deviations from local equilibrium. The former can be expressed as where e and n are the energy density and baryon density in the local rest frame, u µ is the four-velocity of the fluid element (normalized as u 2 = 1), p the pressure given by the EoS, p = p(e, n), and ∆ µν ≡ g µν − u µ u ν . The local rest frame in this work is chosen as the Landau frame specified by the Landau matching conditions [74] which implies u µ n µ = 0 and u µ Π µν = 0. The dissipative term Π µν in the energy-momentum tensor can be written as Π µν = −Π∆ µν + π µν , where Π is the bulk viscous pressure, and π µν the shear stress tensor. The dissipative term n µ in the net baryon current describes its non-zero spatial components in the local rest frame of the fluid. The evolution of these dissipative terms are governed by both microscopic and macroscopic physics, and thus their equations of motion can not be obtained directly from conservation laws. We use the evolution equations from the Denicol-Niemi-Molnar-Rischke (DNMR) theory [75,76], which uses the methods of moments of the Boltzmann equation. In this work, to isolate the effects from baryon diffusion current n µ clearly, we shall ignore the dissipative effects from π µν and Π, focusing only on n µ .
The equation of motion for n µ from the DNMR theory is an Israel-Stewart type equation, whereṅ µ ≡ ∆ µ νṅ ν (the overdot denotes the covariant time derivative D ≡ u µ d µ ), κ n is the baryon diffusion coefficient (also referred to as the baryon conductivity), α ≡ µ/T is the chemical potential µ in the unit of temperature T , and τ n is the relaxation time, on whose scale baryon current relaxes towards its Navier-Stokes limit: here ∇ µ ≡ ∆ µν d ν is the spatial gradient in the local rest frame. The term J µ contains higher order gradient contributions [75,76]. Rewriting Eq. (5) as a relaxation equation,ṅ shows that ∇ µ α is the driving force for baryon diffusion while κ n controls the strength of the baryon diffusion flux arising in response to this force. τ n characterizes the response time scale. We note that both τ n and κ n depend on the microscopic properties of the medium, which have been calculated in various theoretical frameworks, including kinetic theory [52] and holographic models [77,78]. In principle, they can also be constrained phenomenologically by data-driven model inference, but as of today such studies are still very limited for baryon evolution. Rewriting Eq. (5) more specifically for baryon diffusion, we arrive at where θ ≡ d · u, the n µ θ-term arises from J µ (the only term we keep from J µ as given in [54]), δ nn is the associated transport coefficient, and the last two terms come from rewritingṅ µ explicitly, with Γ µ αβ being the Christoffel symbols. Eq. (8) is the equation we use to evolve the baryon diffusion current in this work. We remark that the Navier-Stokes limit of the baryon diffusion current can be rewritten in terms of density and temperature gradients, where the two coefficients are Here χ ≡ (∂n/∂µ) T is the isothermal susceptibility, and w = e + p is the enthalpy density. We note that the gradient expansion is not unique, and writing it in different ways can be used to explore individual contributions separately (see App. C 2). We will discuss the benefits of decomposing the Navier-Stokes limit as Eq. (9) where critical singularities manifest themselves in Sec. III. For later convenience of discussing critical behavior, we also introduce the heat diffusion coefficient, where c p ≡ nT (∂m/∂T ) p is the specific heat, with m ≡ s/n, i.e., the entropy per baryon density [34,79]; λ T is the thermal conductivity, which can be related to the baryon diffusion coefficient κ n by Using Eq. (11) one can also relate the heat diffusion coefficient D p to D B by The above system of hydrodynamic equations is closed by the EoS, either in the form p(e, n) or, equivalently, through the pair of relations µ(e, n), T (e, n) . Again, the EoS is controlled by microscopic physics. We discuss the EoS used in this study in Sec. IV below.

III. CRITICAL BEHAVIOR
In this section we focus on the effects of the QCD critical point on baryon transport in a relativistic QCD fluid. Critical phenomena, albeit ubiquitous, can be classified into certain universality classes determined by the effective degrees of freedom and symmetries of the system. It has been well argued that the QCD critical point belongs to the static universality class of the 3-dimensional Ising model [13,14], and the dynamical universality class of Model H in the Hohenberg-Halperin classification [23,79]. It is also believed that the critical point, if it exists, is beyond the reach of first-principles approaches such as Lattice QCD, and consequently little else is known beyond its universality classification.
Much work has been done to lift the fog from the physics lurking beyond the above universality argument. First, a robust construction of a family of EoS exhibiting the appropriate universal critical properties and matching to existing Lattice QCD calculations has been proposed [65,67]. The EoS encodes the microscopic properties of the QCD matter and is indispensable for solving the macroscopic hydrodynamic equations. Second, a hydrodynamic framework incorporating fluctuations and critical slowing-down has been established, in order to overcome the breakdown of hydrodynamics near the critical point. A deterministic framework, known as hydrokinetic theory, extends conventional hydrodynamics by consistently including fluctuations as additional dynamic degrees of freedom (modes) [80][81][82][83]. The feedback of fluctuating modes renormalizes the bare hydrodynamic variables and gives rise to a delayed response in the form of so-called "long-time tails". The hydro-kinetic approach can be implemented in the critical regime where the fluctuating modes relax to equilibrium on parametrically long time scales [34,40] (see also the reviews [46,84]).
Of course, the inclusion of fluctuations does not by itself cure pathological issues such as acausality or instability of the underlying hydrodynamic framework that arise when straightforwardly extending the nonrelativistic Navier-Stokes equations into the relativistic domain [85,86]. The most-widely used resolution of these issues follows the approach pioneered in Ref. [87], by elevating the dissipative components of the energymomentum tensor to dynamical degrees of freedom subject to their own relaxation-type equations (for which Eq. (7) for the baryon diffusion current is an example). In this approach the causality and stability conditions can be shown to continue to hold in the proximity of the critical point (see App. A).
In the following subsections we discuss how we implement the static and dynamic universal critical behavior, using a simplified setup. Possible future improvements using a more realistic implementation will be discussed in the conclusions (Sec. VI).
A. Implementation of static critical behavior One significant feature of critical phenomena is that, when the system approaches a critical point adiabatically, the equilibrium correlation length, which is typically microscopically small, becomes macroscopically large and eventually diverges. With the purpose of identifying qualitative signatures of a critical point, we characterize all equilibrium quantities exhibiting critical behavior in terms of their parametric dependence on the correlation length. 1 We parametrize the correlation length as follows: Here ξ 0 (µ, T ) is the non-critical correlation length (measured far away from the critical point) while ξ max is an infrared cutoff regulating the divergence at the critical point by implementing a maximum value for the correlation length. The crossover between the critical and non-critical regimes is characterized by the hyperbolic In the above expression (T c , µ c ) is the location of critical point, and ∆µ and ∆T characterize the extent of the critical region along the µ and T axes of the phase diagram; α 1 is the angle between the crossover line (h = 0 axis in the Ising model) and the negative µ axis (see Fig. 1); ν = 2/3, β = 1/3 and δ = 5 approximate the critical exponents of the 3-dimensional Ising universality class [13,14]. Eq. (13) is designed to ensure the following properties: 4) ξ = ξ 0 when |µ − µ c | ∆µ and/or |T − T c | ∆T .
To limit the number of free parameters, we ignore the T and µ dependence of the non-critical correlation length ξ 0 and parametrize the crossover line as [88] T (µ) where T 0 = 155 MeV and κ 2 = 0.0149 are the transition temperature and the curvature of the transition line T (µ) at µ = 0. Location of the critical point (T c , µ c ) is assumed to be on the crossover line [65], and as a consequence Thus T c and α 1 are determined once µ c is provided. Based on the above discussion we choose the following parameter values: parameters, we visualize the correlation length as function of (µ, T ) in Fig. 2.
Several comments are in order: First, our parametrization of the correlation length applies to the crossover region in the left part of the QCD T -µ phase diagram, at µ < µ c , and not to the presumed first-order phase transition at µ > µ c where the theoretical description is complicated by possible phase coexistence and metastability [63]. This suggests choosing the collision beam energy sufficiently high to avoid the latter situation, but not too high to be far from the critical point. Motivated by experimental hints [89] and earlier theoretical studies [52,72] we here set √ s NN = 19.6 GeV. Second, although the correlation length diverges in the thermodynamic limit, heavy-ion collisions create small, rapidly expanding QGP droplets in which finite-size and finitetime effects as well as the critical slowing down [10,27] prevent the correlation length from growing to infinity. A robust estimate for the largest correlation length the system might achieve in this dynamical environment is about 3 fm [27]. The system will thus never get close to our static infrared cutoff ξ max = 10 fm, and our final predictions turn out not to be sensitive to the precise value of this cutoff. Once all thermodynamic quantities and transport coefficients (introduced in the following subsection) are parametrized in terms of ξ as given in Eq. (13), they are defined in both the non-critical and critical regions and thus ready for use in dynamical simulations describing the trajectory of the QGP fireball through the phase diagram. For economy we include in the following discussion not only the dynamic (transport) coefficients but also the thermal susceptibility χ and the specific heat c p which are static (thermodynamic) coefficients.
Near the critical point, fluctuations at the length scale ξ significantly modify the physical transport coefficients, giving rise to their correlation length dependence. In Model H, the shear stress tensor and bulk viscous pressure play important roles in critical dynamics through fluid advection [21,22]. We shall only focus on critical dynamics arising from fluctuations in the hydrodynamic regime (i.e., carrying small wavenumbers/frequencies, q < ∼ 1/ξ), as large-wavenumber fluctuations equilibrate fast compared to the hydrodynamic evolution rate ω ∼ c s k, with c s being the speed of sound. Feedback from off-equilibrium fluctuations that are non-analytic in ω or k, commonly referred to as long-time tails, is suppressed by phase space [21,22] and will be neglected. In other words, the scaling parametrizations below derive from equilibrium fluctuations for thermodynamic quantities and from analytic non-equilibrium fluctuations for transport coefficients.
In the presence of a bulk viscous pressure, the critical contribution to bulk viscosity diverges as ξ z where z = 3 for the QCD critical point [23]. Besides, the associated relaxation time for the bulk viscous pressure also diverges as ξ 3 [72]. In this case, the relaxation rate of the fluctuation modes contributing to bulk viscous pressure is much smaller than the typical hydrodynamic frequency, and they can no longer be treated hydrodynamically, requiring instead an extended framework such as Hydro+ [34,40]. In this work we focus on effects from baryon diffusion and neglect bulk and shear stress effects on the bulk evolution dynamics of the fireball; however, we use the same critical scaling laws for the transport coefficients in Model H as if they were restored (i.e., with non-vanishing shear viscosity η and bulk viscosity ζ).
The following second-order thermodynamic coefficients (isothermal susceptibility χ and specific heat c p ) as well as the first-order transport coefficients (baryon diffusion coefficient κ n and thermal conductivity λ T ) scale with the correlation length as [23] where the exponents are rounded to their nearest integers for simplicity. Therefore, according to Eqs. (9b) and (10), In this work we apply the following parametrizations: where ξ 0 is the non-critical correlation length, κ n,0 is the non-critical value of baryon diffusion coefficient (see Eq. (37) below), χ 0 is the isothermal susceptibility evaluated in the non-critical region, χ 0 ≡ (∂n 0 /∂µ) T where n 0 is the non-critical baryon density. 2 With Eq. (20) the parametrizations with critical scaling for D B , D T and λ T are readily obtained from Eqs. (9b) and (11). We now turn to the critical behavior of the relaxation time τ n . It is worth remembering that the Israel-Stewart type equations (cf. Eq. (5)) provide an ultraviolet completion of the naive (Landau-Lifshitz) hydrodynamic theory. The microscopic relaxation times associated with the new dissipative dynamical degrees of freedom (such as, in our case here, the baryon diffusion current n µ ) play the role of ultraviolet regulators which modify the short-distance (high-frequency) behavior of the theory. For the baryon diffusion current n µ , τ n characterizes the relaxation time to its Navier-Stokes limit n µ NS (which is zero in a homogeneous background). Since n µ can only equilibrate as long as all fluctuating degrees of freedom contributing to n µ also equilibrate, τ n can be considered as the typical equilibration time scale of the slowest fluctuation mode near the critical point. Indeed, in Hydro+/++, the non-hydrodynamic slow-mode evolution equations for critical fluctuations with typical momenta q ∼ ξ −1 have the similar structure as the Israel-Stewart relaxation equations for the dissipative flows arising from thermal fluctuations with wavenumbers ξ −1 ∼ T (see Sec. III C for detailed discussion). Here we assume the scale hierarchy adopted in [40], i.e., the fluctuation wavenumber q is much bigger than the gradient wavenumber k, but still much small than the inverse of thermal length T , i.e., k q T . As already mentioned, Israel-Stewart type equations neglect the non-analytic contributions from long-time tails which we argued above to be negligible (see also Sec. III C).
According to Ref. [40], the slowest mode contributing to n µ is the diffusive-shear correlator between the entropy per baryon density fluctuations δm ≡ δ(s/n) and the flow fluctuations δu µ , i.e., G mµ ∼ δmδu µ . The relaxation rate for this mode with wavenumber q is given by Γ G (q) = (γ η + D p )q 2 where γ η = η/w. The two contributions to this rate stem from the relaxation of the shear stress and of the baryon diffusion, respectively. Near the critical point Γ G is dominated by contributions with typical wavenumbers q ∼ 1/ξ. Given D p ∼ ξ −1 and approximately While the notational distinction between χ and χ 0 is needed here for clarity, we generally drop the subscript "0" for thermodynamic quantities away from the critical region elsewhere to avoid clutter. 3 Another mode contributing to the baryon diffusion current is the pressure-shear mode Gpµ ∼ δpδuµ [40]. Its relaxation rate at wavenumber q is (γ ζ + 4 3 γη + γp)q 2 where γ ζ = ζ/w, γη = η/w, and γp = κnc 2 s T w(∂α/∂p) 2 m . In the presence of bulk viscosity (as we assume in order to ensure the correct scaling in Model H), the relaxation rate for Gpµ is dominated by γ ζ q 2 | q∼ξ −1 ∼ ξ considering γ ζ ∼ ζ ∼ ξ 3 , which is much faster than the rate for the diffusive-shear mode which scales like Γ G ∼ ξ −2 . Even in the absence of the viscosities (i.e. for η = ζ = 0), the dominated Thus it is natural to expect τ n ∼ τ G ∼ ξ 2 . We therefore parametrize τ n as where τ n,0 is the non-critical relaxation time (given explicitly below in Eq. (38)). As we shall discuss in App. A, the parametrization (21) ensures causality.
As an aside, let us comment on the consequences, had we tried to ensure the absence of shear stress by demanding that η = 0. In this case Γ G (q) = D p q 2 | q∼ξ −1 ∼ ξ −3 , and therefore τ n = Γ −1 n ∼ τ G ∼ ξ 3 , which is larger compared to that in the case with shear stress. This arises from the fact that, near the critical point, the shear mode (δu µ ) relaxes to equilibrium parametrically faster than the diffusive mode (δm). As a result, its dissipation changes the scaling exponent of the relaxation time of the diffusive-shear-mode correlator.
Substituting Eq. (13) into Eqs. (20) and (21) we arrive at a complete set of relevant thermodynamic quantities and transport coefficients (i.e., χ, κ n and τ n ) as explicit functions of T and µ that hold in the entire crossover domain of the QCD phase diagram, both far away from and within the critical region.

C. Connection to Hydro+
We close this section by discussing the connection between the Israel-Stewart-like second-order hydrodynamic equations (cf. Eq. (7)) and the evolution equations proposed in the Hydro+ formalism [34].
When comparing the conventional Israel-Stewart formalism without critical effects with the Hydro+ framework, we note that both approaches add nonhydrodynamic degrees of freedom to the conventional hydrodynamic framework, together with their associated relaxation time scales that are not negligible in rapidly evolving systems. The Hydro+ formalism distinguishes itself by the fact that the dynamics of some of these additional non-hydrodynamic modes, the critical slow modes, are controlled by a separate small parameter (different from the conventional Knudsen number for the hydrodynamic gradient expansion) that characterizes the relative "slowness" of the relaxation of the critical slow modes when compared with standard dissipative effects arising from thermal fluctuations.
In our study the critical scaling for the transport coefficients is included in the Israel-Stewart formalism, and consequently it becomes directly comparable to the Hy-dro+ framework with a single wavenumber (q ∼ ξ −1 ) Therefore, the contribution from the pressure-shear mode, which is not the slowest, can be neglected. slow mode. To derive the former from the latter, we introduce a non-hydrodynamic slow mode described by a vector field, in contrast to the scalar field considered as a primary example in the Hydro+ formalism [34]. 4 We denote this field by φ µ and demand that it is a transverse vector, i.e., u · φ = 0, for the sake of later convenience. This vector field will be treated as the slowest mode contributing to n µ , corresponding to G mµ ∼ δmδu µ discussed above.
Including this field the first law of thermodynamics is generalized as follows: Here and below the subscript (+) labels the generalized quantities by taking into account the additional contributions arising from the field φ µ (more precisely, from its deviation from its equilibrium valueφ µ ). The variable π µ is thermodynamically conjugate to φ µ , playing the role of a generalized thermodynamic potential with the constraint u · π = 0. β (+) and α (+) are the associated inverse temperature and chemical potential in units of the temperature, respectively. We require that in thermal equilibrium, where the slow mode φ µ reaches its equilibrium valueφ µ , the entropy density is maximized and equal to its standard equilibrium value In other words, deviations of the slow mode φ µ from its equilibrium valueφ µ reduce the entropy for a given hydrodynamic cell.
The form of F µ φ and A φ shall be specified below. In general, φ µ can also change in response to additional external forces (indicated by · · · in Eq. (24)), such as collective expansion of the background (the θ-term in Eq. (8) arising from J µ ), long-range electromagnetic fields, etc. Here, we focus on its change in response to the gradient of chemical potential in units of the temperature, α (+) , as relevant for our study of baryon diffusion.
The generalized (partial-equilibrium) entropy current is given by where ∆s µ describes a spatial non-equilibrium entropy current in the local rest frame. Using hydrodynamic equations we arrive at where Eqs. (22) and (24) were used. The second law of thermodynamics requires the above expressions to be positive semidefinite, resulting in the following constraints: Here γ π and κ n(+) are (positive semidefinite) transport coefficients. We note that in Eq. (27c) A φ π µ corresponds to the contribution to n µ arising from φ µ . Thus κ n(+) amounts to the baryon transport coefficient in the absence of φ µ (i.e., κ n(+) = κ n,0 ), and n µ approaches the conventional Navier-Stokes limit when φ µ is ignored. The constraint on Π µν (not displayed) leads to its conventional Navier-Stokes form. The last term in (26) can in general have either sign; to always satisfy the second law of thermodynamics we must require it to vanish, giving rise to Eq. (27d). The extended entropy can be decomposed as where we postulate the longitudinal entropy correction due to the mode φ µ as (motivated by, e.g., Ref. [87]) where π φ is the susceptibility. Eq. (29) is quadratic in the space-like vector φ µ and always decreases the entropy. Using the Landau-Khalatnikov formula, the susceptibility π φ can be rewritten as where Γ φ = γ π π φ is the relaxation rate of φ µ , and κ φ is introduced as a new transport coefficient whose physical meaning shall become clear right away. Substituting the above expressions into Eq. (24) and ignoring other external force terms we finḋ Now one can choose φ µ to represent the slowest mode contributing to n µ , i.e., G mµ ∼ δmδu µ and, near the critical point, think of φ µ as the "critical sector" of the baryon diffusion current. As discussed in Sec. III B, for φ µ ∼ G mµ with a typical wave number q ∼ ξ −1 , Γ φ ∼ ξ −2 and κ φ ∼ ξ; near the critical point, it controls the relaxation of the slowest mode contributing to the baryon diffusion. Considering the Navier-Stokes limit of the baryon diffusion given by Eq. (27c), we can write down the relaxation equation for n µ . It receives a contribution of typical wave number q ∼ ξ −1 from φ µ , where we have used the fact that Γ n Γ φ ∼ ξ −2 near the critical point. Here κ n ≡ κ n,0 + κ φ , with κ n,0 denoting the non-critical value of the baryon transport coefficient. Near the critical point κ φ ∼ ξ dominates κ n and consequently κ n κ φ ∼ ξ. The parametrization in Eq. (20) is designed to reproduce this behavior approximately. Finally α (+) is the normalized chemical potential with modifications from φ µ which generate critical behavior in the Equation of State near the critical point. One sees that the single-mode Hydro+ equation (32) (which includes only a single wave number q ∼ ξ −1 ) matches the Israel-Stewart type equation (7) for n µ when the critical scaling (20,21) in the critical regime is accounted for. 5 Eq. (32) can be solved in frequency (ω) space as where is the frequency-dependent baryon diffusion coefficient, with and κ n ≡ κ n (ω=0) being the frequency-independent part. While the imaginary part of Eq. (35) relates the baryonic analog of the electric permittivity, its real part gives rise to the frequency dependence of the baryon transport coefficient: One can infer from Eq. (36) that at small ω, κ n (ω) − κ n (0) ∼ κ n (0) ω 2 while at large ω, κ n (ω) ∼ κ n (0) ω −2 , different from the results given in Ref. [40]. 6 Eqs. (33)-(34) generalize the Navier-Stokes solution for the baryon diffusion current into the critical region where, due to critical slowing down, the slow critical fluctuation modes are not in thermal equilibrium. That is to say, these equations indicate that at large frequencies ω/Γ n > ∼ 1 baryon diffusion is suppressed, hence a naive extrapolation of hydrodynamics with the frequencyindependent coefficient κ n (0) ∼ ξ (see Eq. (18)) to this regime would overestimate the amount of baryon diffusion. Knowing from Eq. (21) that Γ n ∼ ξ −2 , this further implies that the suppression affects modes with frequencies ξ −2 < ∼ ω ξ −1 , i.e. inside the Hydro++ regime discussed in Ref. [40]. When analyzed within the Hydro+ framework with only a single slow mode, the switchingoff of critical contributions to baryon diffusion occurs at ω ξ 2 ∼ 1; this is consistent with the general analysis in Ref. [40] where the full wave number spectrum is taken into account. In other words, the switching-off of the critical contribution to baryon diffusion at finite frequency is taken care of by Eq. (32) (and similarly by Eq. (7)).
Summarizing briefly, we emphasize that in the critical regime the Hydro+ equation (32) leads to similar dynamics as the Israel-Stewart equation (7) for the diffusion current n µ since both equations correctly account for critical slowing down through Γ n = 1/τ n . This correspondence is expected since in the critical regime the off-equilibrium effects are dominated by fluctuations that are effectively frozen. However, Eq. (32) still differs from the Hydro++ equations presented in Ref. [40], since only a single representative mode (with the typical wave number q ∼ ξ −1 ) is analyzed. The exact asymptotic suppression behavior as well as the non-analytic frequency dependence of κ n (cf. footnote 6) arising in Hydro++ from the phase-space integration over all critical modes are both missed by Eq. (32). Indeed, the critical effects on κ n are overestimated by Eq. (32) at small ω (i.e., κ n is less suppressed compared to Ref. [40]). However, we shall see in the following that the resulting overestimation of critical effects is negligible when compared to much stronger suppression effects arising from different origins. For these reasons the single-mode Hydro+ formalism (or, equivalently, the generalized Israel-Stewart formalism amended by critical scaling) serves as a good prototype of the state-of-the-art Hydro+/++ theoryit is sophisticated enough to capture the phenomenologically important feature of critical slowing down while preserving computational economy.

IV. SETUP OF THE FRAMEWORK
In this section we set up the framework for simulating the evolution of a fireball close to the QCD critical point. The core of our framework is the hydrodynamic equations discussed in Sec. II. This requires specification of the EoS and the transport coefficients, as well as initial and final conditions. We also discuss the particlization process of converting the fluid dynamic output into par-ticles whose momentum distributions (after integrating over the conversion hypersurface) can be compared with experimental measurements.
Initial conditions. We start with the initial conditions which, from a physics perspective, describe the initial state of the systems while mathematically providing the initial data for solving the initial value problem associated with our coupled set of partial differential equations. At collision energies of order √ s NN ∼ O(10) GeV, the longitudinal interpenetration dynamics of the two colliding nuclei becomes complicated, in principle requiring a time-dependent, (3+1)-dimensional description of the initial energy-momentum deposition and baryon number doping processes that produce the QGP fluid.
Recently, several so-called "dynamical initialization" algorithms have been proposed to address this problem (see Refs. [47,48,51,90] and the reviews [45,46]). In this exploratory study, however, we try to establish a basic understanding of baryon diffusion dynamics for Au-Au collisions at √ s NN = 19.6 GeV in which we focus entirely on the longitudinal dynamics, modeling a (1+1)dimensional system without transverse gradients initiated instantaneously at a constant proper time τ i (see also Refs. [72,91]). More specifically, we evolve the system hydrodynamically using the longitudinal initial profiles e(τ i , η s ), n(τ i , η s ) provided in Ref. [52], starting at τ i = 1.5 fm/c. 7 The initial hydrodynamic profiles are shown as gray curves in Fig. 4 below. The initial energy density has a plateau covering the space-time rapidity η s ∈ [−3.0, 3.0] whereas the initial net baryon density features a double peak structure and covers a narrower region η s ∈ [−2.0, 2.0], reflecting baryon stopping. 8 For the initial longitudinal momentum flow we take the "static" flow profile u µ = (1, 0, 0, 0) in Milne coordinates (corresponding to Bjorken expansion [92] in Cartesian coordinates), and the initial baryon diffusion current is assumed to vanish, n µ = (0, 0, 0, 0). EoS and transport coefficients. Given these initial conditions, the hydrodynamic equations (1) and (8) are solved by BEShydro [54]. For the EoS at nonzero net baryon density we use neos [53] which was constructed by smoothly joining Lattice QCD data [93][94][95][96] with the hadron resonance gas model. As already mentioned, we here focus on baryon diffusion dynamics by ignoring shear and bulk viscous stresses. For the transport coefficients related to baryon diffusion we rely on 7 Ref. [52] provides longitudinal initial distributions for the entropy and baryon densities. We here adopt the functional form of their initial entropy profile as our energy profile, after appropriate normalization. 8 We note that baryon stopping affects the initial momentum rapidity y of the baryon number carrying degrees of freedom and is typically modelled by a rapidity shift ∆y ∼ 1 − 1.5, depending on system size and collision energy. To translate this rapidity shift into a shift in space-time rapidity ηs (as done in Fig. 4) requires a dynamical initialization model. Different such models yield different initial density and flow profiles [45-48, 51, 90]. 3. Initial longitudinal distributions of κn,0, corresponding to the initial profiles in Fig. 4  the theoretical work in Refs. [52,77,78,97] since phenomenological constraints are still lacking. Specifically, we here use the coefficients obtained from the Boltzmann equation for an almost massless classical gas in the relaxation time approximation (RTA) [52], which gives for the baryon diffusion coefficient

�� ������� �� ���
and for the relaxation time where C n is a free unitless parameter. 9 Throughout this paper, we set C n = 0.4; in Ref. [52] this value was shown to yield good agreement with selected experimental data. Following the kinetic theory approach [52] we also set δ nn,0 = τ n,0 in Eq. (8) as its non-critical value. In the limit of zero net baryon density, κ n,0 remains non-zero at non-zero temperature, lim µ→0 κ n,0 /τ n,0 = nT /(3µ) [52] -a feature also seen in holographic models; for example, using the AdS/CFT correspondence, the (baryon) charge conductivity of r-charged black holes translates into [77,91] κ n,0 = 2π Since κ n,0 is such an important parameter in our study, we offer some intuition about its key characteristics in Fig. 3, where its initial space-time rapidity profile is plotted for three of the theoretical approaches referenced above. 10 The figure shows that in the region |η s | < ∼ 2, where the net baryon density is non-zero (cf. Fig. 4d below), differences exist in κ n,0 (η s ) between the weakly and strongly coupled approaches: In the holographic approaches κ n,0 is suppressed by baryon density while in the kinetic approach it is enhanced. On the other hand we see that at large rapidity |η s | > ∼ 2.5, where the baryon density approaches zero (cf. Fig. 4d), the different models for κ n,0 yield similar distributions, all of them rapidly decreasing towards zero as the temperature decreases (cf. Fig. 4b).
This also implies that generically, as the fireball expands and cools down, all three approaches yield rapidly decreasing amounts of baryon diffusion, as will be discussed below in Sec. V D.
Particlization. After completion of the hydrodynamic evolution (results of which will be discussed in Sec. V) we compute the particle distributions corresponding to the hydrodynamic fields on the freeze-out surface Σ, using the Cooper-Frye formula [99]: where E, p are the energy and momentum of a particle of species i in the observer frame, g i is the spin-isospin degeneracy factor, d 3 Σ µ (x) is the outward-pointing surface normal vector at point x on the three-surface Σ, f eq,i is the equilibrium distribution for particle of species i, and δf diff,i the off-equilibrium correction resulting from net baryon diffusion (other viscous corrections are neglected in this paper). At first order in the Chapman-Enskog expansion of the RTA Boltzmann equation, the dissipative correction from net baryon diffusion δf diff,i is given by [52,100,101] where θ i = 1 (−1) for fermions (bosons), b i is the baryon number of particle species i, p µ ≡ ∆ µν p ν , andκ = κ n /τ n . We will discuss how the critical correction is included, as well as its effects on the final particle distributions, in Sec. V C 2. We evaluate the continuous momentum distribution (40) numerically using the iS3D particlization module [100], ignoring rescattering among the particles and resonance decays after particlization.

V. RESULTS AND DISCUSSION
In this section we discuss the dynamics of the fireball at a fixed collision energy of √ s NN = 19.6 GeV. First we study in Sec. V A baryon diffusion effects on its T -µ trajectories through the phase diagram for cells located at different space-time rapidities η s , and in Sec. V B on freeze-out surface and final particle distributions, in the absence of critical dynamics. Then, in Sec. V C, we discuss how the critical behavior described in Sec. III modifies this dynamics for cells whose trajectories pass close to the critical point. In Sec. V D we point out some generic features of the time evolution of baryon diffusion.

A. Longitudinal dynamics of baryon evolution
In Figure 4, we show snapshots of the longitudinal distributions of the hydrodynamic quantities at two times, the initial time 1.5 fm/c (gray solid lines) and later at τ = 5.5 fm/c, with and without baryon diffusion (red solid and blue dashed lines, respectively). The gray curves in panels (a) and (d) show the initial energy and baryon density distributions from Ref. [52]. The gray lines in panels (b) and (e) show the corresponding temperature and chemical potential profiles, extracted with the neos equation of state [53,[93][94][95][96]. The gray horizontal lines in panels (c) and (f) show the zero initial conditions for the longitudinal flow and baryon diffusion current. The temperature profile (panel (b)) shares the plateau with energy density (panel (a)), up to small structures caused by the double-peak structure of the baryon density (panel (d)) and baryon chemical potential profiles (panel (e)). The chemical potential in panel (e) inherits the double-peak structure from baryon density in panel (d). These structures are also reflected in the pressure (not shown). 11 We next discuss the blue dashed lines in Fig. 4 showing the results of ideal hydrodynamic evolution. Work done by the longitudinal pressure converts thermal energy into collective flow kinetic energy such that the thermal energy density e decreases faster than 1/τ (panel (a)). Small pressure variations along the plateau of the 11 In the very dilute forward and backward rapidity regions one observes a steep rise of the initial µ/T . This feature is sensitive to the rates at which e and n approach zero as |ηs| → ∞, and it is easily affected by numerical inaccuracies. Since both T and µ are close to zero there, the baryon diffusion coefficient κ n,0 also vanishes, and (as seen in Fig. 4f) the apparently large but numerically unstable gradient of µ/T at large ηs does not generate a measurable baryon diffusion current. distribution (caused by the rapidity dependence of µ/T ) lead to slight distortions of the rapidity plateau of the energy density as its magnitude decreases. Longitudinal pressure gradients at the forward and backward edges of the initial rapidity plateau accelerate the fluid longitudinally, generating a non-zero η s -component of the hydrodynamic flow at large rapidities (panel (c)). As seen in panels (a) and (c), the resulting longitudinal rarefaction wave travels inward slowly, leaving the initial Bjorken flow profile u η = 0 untouched for |η s | < 2.5 up to τ = 5.5 fm/c. For Bjorken flow without transverse dynamics baryon number conservation implies that nτ remains constant. Panel (d) shows this to be the case up to τ = 5.5 fm/c because, up to that time, the initial Bjorken flow has not yet been affected by longitudinal acceleration over the entire η s -interval in which the net baryon density n is non-zero. Panel (e) shows, however, that in spite of nτ remaining constant within that η s range, the baryon chemical potential µ/T decreases with time, as required by the neos equation of state.
The nontrivial evolution effects of turning on the baryon diffusion current via Eq. (8) are shown by the red solid lines in Fig. 4. The baryon diffusion current itself is plotted in panel (f) and will be discussed shortly. Panels (a)-(c) show that baryon diffusion has almost no effect at all on the energy density (and, by implication, on the pressure), the temperature, and the hydrodynamic flow generated by the pressure gradients. Given the weak dependence of pressure and temperature on baryon density through the EoS this is to be expected. Baryon diffusion does, however, significantly modify the rapidity profiles of the net baryon density (d), chemical potential µ/T (e). Generated by the negative gradient of µ/T , the baryon diffusion current moves baryon number from high-to low-density regions, causing an overall broadening of the baryon density rapidity profile in (d) while simultaneously filling in the dip at midrapidity [48,51,52,54,56,91]. Panel (e) shows how the chemical potential µ/T tracks these changes in the baryon density profile (panel (d)), and panel (f) shows the baryon diffusion current responsible for this transport of baryon density, with its alternating sign and magnitude tracing the sign and magnitude changes of −∇(µ/T ). As we shall see in Sec. V D, the smoothing of the gradients of baryon density and chemical potential contributes to a fast decay of baryon diffusion effects. Fig. 4 indicates non-trivial thermal, chemical, and mechanical evolution at different rapidities. Fluid cells at different η s pass through different regions of the QCD phase diagram and may therefore be affected differently by the QCD critical point [72,102,103]. This has led to the suggestion [104] of using rapidity-binned cumulants of the final net proton multiplicity distributions as the possibly sensitive observables of the critical point. 12 To illustrate the point we show in Fig. 5 the phase diagram 12 We caution that at BES energies the mapping between space- trajectories of fluid cells at several selected |η s | values, 13 both with and without baryon diffusion. As we move from mid-rapidity to |η s | = 2.0, the starting point of these trajectories first moves from µ 0.28 GeV at η s = 0 to the larger value µ 0.45 GeV at η s = 1.5, but then turns back to µ 0.2 GeV at η s = 1.75, and finally to µ 0 at η s = 2.0, without much variation of the initial temperature T i 0.25 GeV (see Figs. 4b,e). The difference between the dashed (ideal) and solid (diffusive) trajectories exhibits a remarkable dependence on η s : Both the sign and the magnitude of the diffusion-induced shift in baryon chemical potential depend strongly on space-time rapidity. In most cases, we note that the solid (diffusive) trajectories move initially rapidly away from the corresponding ideal ones, but then quickly settle on a roughly parallel ideal trajectory. A glaring exception is the trajectory of the cell at η s = 1.5, which starts at the maximal initial baryon chemical potential and keeps moving away from its initial ideal T -µ trajectory for a long period, settling on a new ideal trajectory only shortly before it reaches the hadronization phase transition. The reason for this behavior can be found in Fig. 4e, which shows that at η s = 1.5 the gradient of µ/T remains large throughout the fireball evolution. But almost everywhere else baryon diffusion effects die out quickly.
time rapidity ηs of the fluid cells and rapidity y of the emitted hadrons is highly nontrivial and requires dynamical modelling. 13 Cells at opposite but equal space-time rapidities are equivalent because of ηs → −ηs reflection symmetry in this work. Since ideal fluid dynamics conserves both baryon number and entropy, the dashed trajectories are lines of constant entropy per baryon. This is shown by the dashed lines in Fig. 6. Baryon diffusion leads to a net baryon current in the local momentum rest frame and thereby changes the baryon number per unit entropy. This is illustrated by the solid lines in Fig. 6. Depending on the direction of the µ/T gradients, baryon diffusion can increase or decrease the entropy per baryon.
We close this discussion by commenting on the turning of the dashed m ≡ s/n = const. trajectories in Fig. 5 from initially pointing towards the lower left to later pointing towards the lower right. This is a well known feature of isentropic expansion trajectories in the QCD phase diagram [53,105,106] that reflects the change in the underlying degrees of freedom, from quarks and gluons to a hadron resonance gas, at the point of hadronization as embedded in the construction of the EoS. Figure 5 is reminiscent of the QCD phase diagram often shown to motivate the study of heavy-ion collisions at different collision energies in order to explore QCD matter at different baryon doping (see, for example, the 2015 DOE-NSF NSAC Long Range Plan for Nuclear Physics [1]). What had been shown there are (isentropic) expansion trajectories for matter created at midrapidity in heavy-ion collisions with different beam energies, whereas Fig. 5 shows similar expansion trajectories for different parts of the fireball in a collision with a fixed beam energy. Fig. 5 thus makes the point that in general the matter created in heavy ion collisions can never be characterized a single fixed value of µ/T . At high collision energies space-time and momentum rapidities are tightly correlated, η s y, and different η s regions with different baryon doping µ/T can thus be more or less separated in experiment by binning the data in momentum rapidity y. This motivates the strategy of scanning the changing baryonic composition in the T -µ diagram by performing a rapidity scan at fixed collision energy rather than a beam energy scan at fixed rapidity [72,102,104]. This strategy fails, however, at lower collision energies where particles of fixed momentum rapidity can be emitted from essentially every part of the fireball and thus receive contributions from regions with wildly different chemical compositions, with non-monotonic rapidity dependences that are non-trivially and non-monotonically affected by baryon diffusion.
B. Freeze-out surface and final particle distributions The expansion trajectories shown in the previous subsection all end at the same constant proper time (see Fig. 6). In phenomenological applications it is usually assumed that the hydrodynamic stage ends and the fluid falls apart into particles when all fluid cells reach a certain "freeze-out energy density", here taken as e f = 0.3 GeV/fm 3 . 14 With such a freeze-out criterion, fluid cells at different η s freeze out at different times τ f (η s ). In this subsection we discuss this freeze-out surface and the distributions of particles emitted from it. Fig. 7 shows the freeze-out surface τ f (η s ) in panel (a) as well as the longitudinal flow, baryon chemical potential, and longitudinal component of the baryon diffusion current in panels (b)-(d). 15 Ideal and diffusive hydrodynamics are distinguished by blue dashed and red solid lines. Panel (a) shows that initially the longitudinal pressure gradient causes the fluid to grow in η s direction before it starts to shrink after τ > ∼ 4 fm/c due to cooling and surface evaporation. As seen in Fig. 4a, the core of the fireball remains approximately boost invariant while cooling by performing longitudinal work, until the longitudinal rarefaction wave reaches it. Once the energy density in this boost-invariant core drops below e f , it freezes out simultaneously, as seen in the flat top of the freezeout surface shown in panel (a). Slight deviations from boost invariance are caused by the effects of the boostnon-invariant net baryon density profile and its (small) effect on the pressure whose gradient drives the hydrodynamic expansion. Baryon diffusion has practically no effect on the freeze-out surface, nor on the longitudinal flow along this surface shown in panel (b), owing to the weak dependence of the EoS on baryon doping. The distributions of the baryon chemical potential and baryon diffusion current across this surface, on the other hand, are significantly affected by baryon diffusion, as seen in panels (c) and (d). It bears pointing out, however, that Given these quantities on the freeze-out surface, we use the iS3D module [100] to evaluate the Cooper-Frye integral (40,41) for the rapidity distributions of hadrons emitted from the freeze-out surface. Results are shown in Fig. 8. Panel (b) indicates that baryon diffusion has negligible effects on meson distributions. It affects only baryon distributions. Panel (a) shows that baryon diffusion significantly increases the proton and net-proton yields at mid-rapidity and also broadens their rapidity distributions at large rapidity. Both of these effects on the baryon distributions were also observed in earlier work that, different from our simulations, additionally included resonance decays and full hadronic rescattering [48,51,52,91]; furthermore, they were found to increase with the magnitude of the baryon diffusion coefficient κ n . The approximate boost-invariance of the longitudinal flow over a wide range of η s on the freeze-out surface (see Fig. 7b) maps the baryon diffusion effects seen in Figs. 4d,e and 7c as functions of space-time rapidity η s onto momentum rapidity y p in Fig. 8a [48,51,52,91]. We take advantage of the iS3D option to include both the Chapman-Enskog and 14-moment approximations for the dissipative corrections (41), comparing the two in Fig. 8. The difference is seen to be negligibly small, and even ignoring in Eq. (40) δf diff,i entirely does not make much of a difference (not shown). This reflects the tiny magnitude of the baryon diffusion current on the freeze-out surface seen in Fig. 7d. 16 We emphasize that the mapping of baryon diffusion effects seen as a function of spacetime rapidity η s in Figs. 4d,e and 7c onto momentum rapidity y p is expected to be model dependent, and may not work for initial conditions in which the initial velocity profile is not boostinvariant or the initial η s -distribution of the net baryon density looks different. This initial-state modeling uncertainty has so far prohibited a meaningful extraction of the baryon diffusion coefficient from experimental data (see, however, Ref. [52] for a valiant effort). Additional uncertainties from possible critical effects associated with QCD critical point on the bulk dynamics, especially through baryon diffusion, may further complicate the picture, in particular as long as the location of the critical point is still unknown. In the following subsection we address some of these effects arising from critical dynamics.

C. Critical effects on baryon diffusion
In this section, we explore whether the QCD critical point can have significant effects on the bulk dynamics, through the baryon diffusion current. For this purpose, we include critical effects as described in Sec. III, and explore effects from critical slowing down on the hydrodynamic transport (Sec. V C 1), as well as critical corrections to final particle distributions through the Cooper-Frye formula (Sec. V C 2).

Critical slowing down of baryon transport
As discussed in Sec. III, in the critical region baryon transport is affected by critical slowing down [23]. Outside the critical region all thermodynamic and transport properties approach their non-critical baseline described in Sec. V A, but as the system approaches the critical point its dynamics is affected by critical modifications of the transport coefficients involving various powers of ξ/ξ 0 > 1. We study this by incorporating the critical scaling of χ, κ n and τ n in Eqs. (20) and (21), with the correlation length ξ(µ, T ) parametrized by Eq. (13).
Before doing any simulations we briefly discuss qualitative expectations. Eq. (19) indicates that, as the correlation length grows, ξ/ξ 0 > 1, the coefficient D B is suppressed while D T is enhanced. According to Eqs. (9) a suppression of D B reduces the contribution from baryon density inhomogeneities while an enhancement of D T increases the contribution from temperature inhomogeneities to the Navier-Stokes limit n ν NS = κ n ∇ ν (µ/T ). 17 In addition to thus moving its Navier-Stokes target value, proximity of the critical point also increases the time τ n (see Eq. (21)) over which the baryon diffusion current relaxes to its Navier-Stokes limit -its response to the driving force is critically slowed down.
Repeating the simulations with the same setup as in Sec. V A, except for the inclusion of critical scaling, yields the results shown in Fig. 9. For the parametrization of the correlation length ξ(µ, T ) we assumed a critical point located at (T c = 149 MeV, µ c = 250 MeV). This is very close to the right-most trajectory shown in Fig. 9 17 We note that in the literature sometimes only the baryon density gradient term D B ∇n is included in the diffusion current (see, e.g., Refs. [23,79]) which then leads to its generic suppression close to the critical point. which should therefore be most strongly affected by it. 18 Surprisingly, none of the trajectories, not even the one passing the critical point in close proximity, are visibly affected by critical scaling of transport coefficients.
To better understand this we plot in Fig. 10 the history of the correlation length and baryon diffusion current at different η s . In panel (a) we see that ξ does show the expected critical enhancement, by up to a factor ∼ 4.5 at η s = 1. This maximal enhancement corresponds to τ n 20 τ n,0 and D B 0.22 D B,0 , naively suggesting significant effects on the dynamical evolution. However, the critical enhancement of the correlation length does not begin in earnest before the fireball has cooled down to a low temperature T < ∼ T c + ∆T . Fig. 10b shows at at this late time the baryon diffusion current has already decayed to a tiny value. 19 In other words, the largest baryon diffusion currents are created at early times when the temporal gradients are highest but the system is far from the critical point; by the time the system gets close to the critical point, thermal and chemical gradients have decayed to such an extent that even a critical enhancement of the correlation length by a factor 5 can no longer revive the baryon diffusion current to a noticeable level. 18 Since we do not have the tools here to handle passage through a first-order phase transition, we do not consider any expansion trajectories cutting the first-order transition line to the right of QCD critical point. 19 This statement remains true if one multiplies n η with the metric factor √ −g = τ to obtain the baryon diffusion current in physical units of fm −3 . This two-stage feature, with a first stage characterized by large baryon diffusion effects without critical modifications and a second stage characterized by large critical fluctuations [21,22] with negligible baryon diffusion effects on the bulk evolution, is an important observation. For a deeper understanding we devote Sec. V D to a more systematic investigation of the time evolution of the diffusion current, but not before a brief exploration in the following subsection of critical effects on the final singleparticle distributions.

Critical corrections to final single-particle distributions
How to consistently include critical fluctuation effects on the finally emitted single-particle distributions, at the ensemble-averaged level, is still a subject of active research. A solid framework may require a microscopic picture involving interactions between the underlying degrees of freedom and the fluctuating critical modes during hadronization [109]. In this work we employ a simple ansatz where critical corrections to the final particle distributions are included only via the diffusive correction δf diff from Eq. (41) appearing in the Cooper-Frye formula (40). In this subsection this dissipative correction is computed from the simulations described in the preceding subsection which include critical correlation effects through critically modified transport coefficients, specifically a normalized baryon diffusion coefficientκ ≡ κ n /τ n with critical scalingκ obtained from Eqs. (20) and (21) usingκ 0 = κ n,0 /τ n,0 . 20 Since we saw in the preceding subsection that the hydrodynamic quantities on the particlization surface are hardly affected by the inclusion of critical scaling effects during the preceding dynamical evolution, the main critical scaling effects on the emitted particle spectra arise from any critical modification thatκ might experience on the particlization surface. The space-time rapidity distribution of the correlation length ξ along the freeze-out surface, as well as the net proton rapidity distributions with and without critical scaling effects, are shown in Fig. 11. Panel (a) shows that ξ peaks near |η s | 1.0 on the freeze-out surface, consistent with Fig. 10a. Note that, although fluid cells at different η s generally freeze out at different times, the freeze-out surface in Fig. 7a shows that within η s ∈ [−1.5, 1.5] all fluid cells freeze out at basically the same time τ f ∼ 17 fm/c. Therefore Fig. 11a indeed corresponds to the ξ values at different η s at the end of the evolution in Fig. 10a.
Even though Fig. 11a shows a critical enhancement of ξ/ξ 0 2.7 near η s 1.0, corresponding toκ/κ 0 0.37, we see in Fig. 11b that the net proton distribution is modified by at most a few percent. The lower panel in Fig. 11b indicates that the largest critical corrections indeed correspond to regions of large ξ/ξ 0 , sign-modulated by the direction of the baryon diffusion current (cf. Fig. 7d). We also notice a thermal smearing when mapping the distribution of ξ in η s to the modification of net proton distribution in y p . The critical modification of the net proton spectra arising from the diffusive correction to the distribution function is very small also because δf diff in (41) is roughly proportional to the magnitude of n µ which is tiny. Such small modifications are certainly unresolvable with current or expected future measurement. Indeed, the magnitudes of these modifications depend on the correlation length on the freeze-out surface which in turn depends on the choice of the freeze-out energy density e f . Independence of the final results from this choice could be achieved by properly sampling the critical correlations on this surface and then propagating them to the completion of kinetic freeze-out with a hadronic transport code that appropriately accounts for critical dynamics; unfortunately, these options are presently not yet available.
In conclusion, critical scaling effects on both the hydrodynamic evolution of the bulk medium and the finally emitted single-particle momentum distributions are small, mostly because by the time the system passes the critical point and freezes out the baryon diffusion current has decayed to negligible levels.

D. Time evolution of baryon diffusion
In this subsection we further analyze the baryon diffusion dynamics and the origins of its rapid decay. We define Knudsen and inverse Reynolds numbers for baryon diffusion and display their space-time dynamics. The resulting insights are relevant for model building and for the future quantitative calibration of the bulk fireball dynamics at non-zero chemical potential.
FIG. 12. Same as Fig. 10b, but zoomed in onto the early evolution stage and including for comparison as dashed lines the corresponding Navier-Stokes limit n η NS of the longitudinal diffusion current.

Fast decay of baryon diffusion
As discussed in Sec. II, the diffusion current relaxes to its Navier-Stokes limit n ν NS = κ n ∇ ν (µ/T ) on a time scale given by τ n . General features of baryon diffusion evolution can thus be understood by following the time evolution of n µ NS , κ n and τ n . Here we focus on their evolution without inclusion of critical scaling since we established that the latter has negligible effect on the bulk evolution and therefore the non-critical values of n µ NS , κ n and τ n evolve almost identically with and without inclusion of critical effects. Fig. 12 shows a comparison of the longitudinal baryon diffusion current (solid lines) with its Navier-Stokes limit (dashed lines) at different space-time rapidities. One sees that the relaxation equation for n η tries to align the diffusion current with its Navier-Stokes value (which is controlled by the longitudinal gradient ∇ η (µ/T )) but the finite relaxation time delays the response, causing n η to perform damped oscillations around n η NS . This is most clearly illustrated in Fig. 12 by following the cell located at η s = 1.5 (uppermost): Initialized at zero, n η initially rises steeply, trying to adjust to its positive and rapidly increasing Navier-Stokes value, but at τ 1.7 fm/c the longitudinal gradient of µ/T switches sign and n η NS starts to decrease again. The hydrodynamically evolving n η follows suit, turning downward with a delay of about 0.3 fm/c (which, according to Fig. 13b below, is the approximate value of the relaxation time τ n,0 at τ = 2 fm/c), but soon finds itself overshooting its Navier-Stokes value. For the cell located at η s = 1.25, n η crosses its Navier-Stokes value even twice.
As long as the relaxation time τ n is short and not dramatically increased by critical slowing down, the rapid decrease of the dynamically evolving diffusion current is seen to be a generic consequence of a corresponding rapid decrease of its Navier-Stokes value: Fig. 12 shows that after τ ∼ 3.5 fm/c, n η basically agrees with its Navier-Stokes limit n η NS . Fig. 13b shows that in the absence of critical effects the relaxation time τ n,0 = C n /T increases by less than a factor of 2 over the entire fireball lifetime. The rapid decrease of n η NS is a consequence of two factors: (i) the gradients of µ/T decrease with time, owing to both the overall expansion of the system and the diffusive transport of baryon charge from dense to dilute regions of net baryon density, and (ii) the baryon diffusion coefficient κ n,0 decreases dramatically (by almost an order of magnitude over the lifetime of the fireball as seen in Fig. 13a), as a result of the fireball's decreasing temperature.
In summary, three factors contribute to the negligible influence of the QCD critical point on baryon diffusion: First, baryon diffusion is largest at very early times when its relaxation time is shortest and it quickly relaxes to its Navier-Stokes value; the latter decays quickly, due to decreasing chemical gradients and a rapidly decreasing baryon diffusion coefficient. Second, the relaxation time for baryon diffusion increases at late times, generically as a result of cooling but possibly further enhanced by critical slowing down if the system passes close to the critical point. This makes it difficult for the baryon diffusion current to grow again. Third, critical effects that would modify 21 the Navier-Stokes limit for the baryon diffusion current become effective only at very late times when n η NS has already decayed to non-detectable levels. The baryon diffusion current thus remains small even if its Navier-Stokes value were significantly enhanced by critical scaling effects.

Knudsen and inverse Reynolds numbers
We close this section by investigating the (critical) Knudsen and inverse Reynolds numbers associated with baryon diffusion. These are typically taken as quantitative measures to assess the applicability of second order viscous hydrodynamics such as the BEShydro framework employed in this work. Copying their standard definitions for shear and bulk viscous effects [52,54,111], we here set for baryon diffusion, where θ is the scalar expansion rate. Kn is the ratio between time scales for microscopic diffusive relaxation (τ n ) and macrosopic expansion (τ exp = 1/θ); the relaxation time τ n includes the effects of critical slowing down in the neighborhood of the QCD critical point. Re −1 is the ratio between the magnitude of the off-equilibrium baryon diffusion current and the equilibrium net baryon density in ideal fluid dynamics. Their space-time evolutions are shown in Fig. 14, together with the freeze-out (particlization) surface at e f = 0.3 GeV/fm 3 . Fig. 14a tells us that Kn > ∼ 1 happens only outside the freeze-out surface, in the fireball's corona where the fluid has already broken up into particles even at the earliest stage of the expansion. The short-lived peak in Kn near τ = τ i = 1.5 fm/c and η s ∼ 4 is caused by the rapid increase of τ n in the dilute and very cold corona of the fireball (note that τ n,0 = C n /T ). Critical slowing down near the QCD critical point causes the Knudsen number to increase somewhat around η s = 1 close to the freeze-out surface; this critical enhancement is barely visible as a light cloud on a blue background, indicating critical Knudsen numbers in the range Kn ∼ 0.5 − 0.7. Fig. 14b, on the other hand, indicates that Re −1 < ∼ 0.3 during the entire evolution, even close to the places where the Navier-Stokes value of the baryon diffusion current peaks at early times (see Fig. 12). After τ ∼ 5 fm/c (including the entire critical region around the QCD critical point) its maximum value drops below 0.1, reflecting of the rapid decay of the baryon diffusion current. The maximal Re −1 occurs around η s 2 shortly after the hydrodynamic evolution starts at τ i = 1.5 fm/c. From it emerges a region of sizeable inverse Reynolds number which ends at two moving boundaries where Re −1 = 0 (dark blue). The left boundary, moving towards smaller η s values, reflects a sign change of the baryon diffusion current (see Fig. 4f where at τ = 5.5 fm/c n η flips sign at η s 1). The right boundary, on the other hand, corresponds to where n µ decays to zero (which, according to Fig. 4f, happens at η s 2.5 when τ = 5.5 fm/c). The initial outward movement of the right boundary is a result of baryon transport to larger space-time rapidity. It stops moving after the diffusion current has decayed and no longer transports any baryon charge longitudinally.
The small values of Kn and Re −1 during the entire fluid dynamical evolution validate the applicability of second order viscous hydrodynamics, BEShydro, for describing flow and diffusive transport of baryon charge in the collision system studied here.

VI. CONCLUSIONS AND DISCUSSION
Baryon diffusion is an important dissipative effect in the hydrodynamic evolution of systems carrying a conserved net baryon charge. It smoothes out chemical inhomogeneities by transporting baryon charge relative to the momentum flow from regions of large to smaller net baryon density. It is driven by gradients of µ/T , but the corresponding transport coefficient characterizing the magnitude of the diffusive response, the baryon diffusion coefficient κ n (µ, T ), is still poorly constrained theoretically. In this work we studied a (1+1)-dimensional system without transverse gradients and flow to explore diffusive baryon transport along the longitudinal (beam) direction in heavy-ion collisions and thereby gain intuition about possible strategies to extract the baryon diffusion coefficient from experimental data. A fundamental difficulty for such an extraction is that baryon diffusion manifests itself as a transport of net baryon density in coordinate space while the experimental observations provide a snapshot of the hydrodynamic medium at the end of its lifetime in momentum space. In a rapidly expanding system, collective flow gradients map different space-time regions into different regions of momentum space, but the thermal momentum spread in the local rest frame blurs this map (less so for the heavier baryons than for the more abundant lighter mesons), and makes it difficult to reconstruct the movement generated by baryon diffusion in coordinate space from the final net baryon distribution in momentum space. In longitudinal direction local thermal motion results in a rapidity spread of order T / m T [112,113] where T is the kinetic freeze-out temperature and m T is the average transverse mass of the particle species in question. The flow-induced separation of different fireball regions in longitudinal momentum rapidity is easier at higher collision energies where the fluid medium created in the collision covers a wider rapidity range, i.e. a larger multiple of the thermal smearing width. At lower beam energies, such as those probed in the RHIC BES campaign, unfolding the measured final rapidity distribution into different regions of space-time rapidity, with possibly different net baryon densities, becomes much harder. And therefore the reconstruction of baryon-diffusion induced matter transport across space-time rapidity also becomes much harder.
In this work we studied central Au-Au collisions at √ s NN ∼ 20 GeV in which the fireball covers about 7 units of space-time rapidity along the beam direction (Fig. 4a). Assuming that baryon stopping leads to a space-time rapidity shift of about 1.5 units for the incoming projectile and target baryons, the initial net baryon distribution had a width of about 4 units. It was modeled by a double-humped function with two well-separated peaks located at η s = ± 1.5 (Fig. 4d). After accounting for ideal hydrodynamic evolution and thermal smearing this resulted in a double-humped net proton rapidity distribution whose peaks were still relatively cleanly separated by about 2 units of rapidity, but after including baryon diffusion they almost (though not quite) merged into a single broad hump around midrapidity (Fig. 8a). In our calculation, the QCD critical point was positioned at a baryon chemical potential µ = 250 MeV. Recent Lattice QCD results put the likely location of this critical point at µ > 400 MeV [96,114,115], which requires lower collision energies for its experimental exploration. At lower collision energies the width of the initial space-time rapidity interval of non-zero net baryon density will be narrower, and the final net-proton distribution will eventually become single-peaked in central collisions [116].
In the work presented here we focused on the questions how diffusive baryon transport manifests itself along the beam direction in hydrodynamic simulations, what traces it leaves in the finally measured rapidity distributions, and how it is affected by critical scaling of transport co-efficients in the proximity of the QCD critical point. To address these questions we systematically discussed the static and dynamic critical behavior of thermodynamic properties (especially those associated with baryon transport) and introduced an analytical parametrization of the correlation length that correctly reproduces the critical exponents of the 3D Ising universality class. Based on a careful comparison with the Hydro+/++ framework [34,40,72] we identified the critical scaling ("critical slowing down") of the relaxation time for the baryon diffusion current (τ n ∼ ξ 2 ), and demonstrated that in the critical regime the Israel-Stewart type equation for the baryon diffusion current plays the role of a single-mode Hydro+ equation for a vector slow mode. We did not discuss the out-of-equilibrium evolution of the slow mode itself. For a single scalar slow mode, this was studied elsewhere using the BEShydro+ framework [22] where it was found that its feedback to the hydrodynamic bulk evolution was negligible [21,22]. A systematic Hydro++ study incorporating the full coupled evolution of all relevant non-hydrodynamic critical slow modes is still outstanding.
We are not the first to point out that, because of its extended nature, different regions within a collision fireball probe different regions of the QCD phase diagram as they cool and expand. We found that, at early times, strong longitudinal gradients of µ/T lead to significant longitudinal baryon diffusion currents that shift the expansion trajectories for different parts of the fireball in different directions within the QCD phase diagram. The final net proton rapidity distributions reflect these shifts, albeit blurred by thermal smearing. Our model assumes a boost-invariant initial longitudinal momentum flow y flow = η s which maximizes the correlation between the space-time rapidity of a fluid cell on the freeze-out surface and the momentum rapidity of the final hadrons it emits. The expected breaking of boostinvariance of the longitudinal flow pattern by dynamical initialization effects at lower collision energies even near midrapidity [47,48,117] will inevitably weaken this correlation, adding to the decorrelating effects of thermal smearing. This will likely result in significant sensitivities of baryon diffusion coefficients inferred from experimental net baryon rapidity distributions to poorly controlled model ambiguities in the initial space-time rapidity profile of the net baryon density assumed in the dynamical model.
The baryon diffusion flows observed in the calculations presented in this paper are characterized by an important feature: They show almost no sensitivity to critical effects even for cells passing close to the critical point. Taken at face value, this implies that the hydrodynamic evolution of baryon diffusion leading to the finally emitted ensemble-averaged single-particle momentum spectra does not carry useful information for locating the QCD critical point.
The absence of critical effects on baryon diffusion in this work contrasts starkly with the strong critical ef-fects on the evolution of the bulk viscous pressure found in Ref. [72] which led to significant distortions of the rapidity distributions for all emitted hadron species. The main reason for that behavior is that the bulk viscosity and bulk viscous pressure are generically large around the quark-hadron phase transition, even without explicitly including critical scaling effects. Adding the latter thus causes significant modifications of the dynamics of the bulk viscous pressure. 22 Baryon diffusion effects, on the other hand, here appear to be insensitivity to critical dynamics, for two main reasons: (i) The baryon diffusion flows are strong at early times but decay very quickly, before the system enters the critical region, because diffusion reduces the initially strong chemical gradients ∇(µ/T ) that drive the flows, and the baryon diffusion coefficient describing the response to these gradients decreases quickly as the fireball cools by expansion. (ii) By the time the system reaches the phase transition, possibly passing close to the critical point, the Navier-Stokes value of the baryon diffusion current is already very small; critical enhancement by the baryon diffusion coefficient κ n by a power of the correlation length ξ/ξ 0 therefore does not help to revive it, and in any case the relaxation rate controlling the approach of the baryon diffusion current to its critically affected Navier-Stokes value is reduced by critical slowing down. The smallness of the baryon diffusion current on the freeze-out surface also implies very small dissipative corrections to the Cooper-Frye formula at particlization.
The observed insignificance of critical effects on baryon diffusion might be taken as permission to calibrate the fireball medium's bulk evolution at BES energies without worrying about the QCD critical point and its location. This may be premature, however, for multiple reasons: (1) The inclusion of shear and bulk viscous effects will modify the expansion trajectories through the phase diagram. Although the critical enhancement of the shear viscosity is negligible (η ∼ ξ /19 [23] where = 4−d, with d being the number of spatial dimensions), critical slowing down of the shear stress tensor may still be significant, possibly causing larger residual shear stresses in the vicinity of the critical point than predicted in the absence of critical effects. This has not been studied. The bulk viscous pressure is already known to be strongly affected by critical phenomena [72]. We see a potential for these effects to spill over into the baryon diffusion channel through second-order transport coefficients that couple these channels. Should this happen, shear and bulk viscous effects could invalidate some of the findings of the work presented here (in particular the rapid decay of the baryon diffusion current observed in Sec. V D 1). To clarify this, a future full simulation should include all dissipative effects simultaneously. (2) Large uncertain-ties still exist about the size of the critical region which is determined by the parametrization of the correlation length [65,67]. (3) At lower beam energies, the system may enter the critical region earlier, before the baryon diffusion significantly decays. When the evolution is fully (3+1)-dimensional the edge of the fireball may enter the critical region earlier as well, while the diffusion current is still appreciable.
A fully quantitative evaluation of the significance of critical effects on bulk medium evolution at BES energies can only be made on top of an at least tentatively constrained bulk evolution that includes all necessary theoretical ingredients and complications [57]. Only if this confirms negligible critical effects on bulk evolution can the final model calibration be safely made without worrying about critical modifications. convention of [65,119,120], the linear mapping from 3dimensional Ising variables (i.e., reduced Ising temperature r and magnetic field h) to the coordinate variables of QCD phase diagram (i.e., temperature T and baryon chemical potential µ) reads where (w, ρ) are scale factors for the Ising variables r and h. For notation simplicity we let s i = sin α i , c i = cos α i , i = 1, 2, where α 1 and α 2 are the angles relative to the negative µ axis for the mapped r and h axis respectively, defined in Fig. 1. From Eq. (B1) one can also find the inverse mapping relations Applying the flow with zero transverse components, these equations can be written into coupled partial differential equations, which can be solved using other softwares (e.g. Mathematica) for testing BEShydro's performance. While these equations clearly exhibit the physics in the LRF (which varies from point to point), BEShydro solves the conservation laws (Eqs. (1a,1b)) in a fixed global computational frame; thus this would be a non-trivial test of our hydrodynamic code.
To solve the equations, we start the system at initial time τ i = 1 fm, using the following initial conditions: u µ = (1, 0, 0, 0) and the initial diffusion current is zero; the longitudinal distribution of baryon density is given as where the initial temperature T (τ i , η s ) = T 0 = 0.5 GeV, and n max = 10 fm −3 , σ η = 0.5 and η ± = ± 1. Here g s = 16 is the degeneracy factor in the EoS below. The initial chemical potential is set to be zero everywhere. One also needs the two transport coefficients in Eq. (C4), and we use [76] κ n = 3 16σ tot , τ n = 9 4 1 nσ tot , with σ tot = 10 mb = 1 fm 2 . To close the equations, we use the EoS e = 3p = 3nT , n = g s π 2 T 3 exp µ T . (C7) We note that the EoS can be analytically inverted to get T (e, n) and µ(e, n), directly usable for hydrodynamic codes, which is convenient for testing such codes. The setup above which is from Ref. [56] can be solved semi-analytically using e.g. Mathematica and used for validating BEShydro, the comparison of which is shown in Fig. 15, for the baryonic sector. One can see from the figure that BEShydro works perfectly well in the longitudinal evolution.
With the solutions, we can also get the space-time distribution of the freeze-out surface, and correspondingly the other hydrodynamic quantities on the surface; thus we can take the chance to test the freeze-out finder implemented in BEShydro, which is based on Cornelius [107]. We have tested its efficiency within BEShydro at non-zero baryon density in the transverse plane [108]. Here, we validate it for this work, where longitudinal dynamics with baryon diffusion current is evolved. In Fig. 16, we define the freeze-out surface at two different freeze-out energies, and we see that the space-time profile and distribution of the freeze-out temperature perfectly agree with the semi-analytical results.

Calculation of various gradients
As mentioned in Sec. II, to include critical contributions to the EoS without using such an EoS, we rewrite the Navier-Stokes limit term of baryon diffusion into two terms: Note that consistence between results from κn∇α (red solid line) and DT ∇T + DB∇n (green dashed line) is a highly nontrivial test of numeric methods, since the later involves calculating (∂p/∂T )n and interpolating tabulated χ etc. Though the overall agreement is very good, small wiggles can be seen near the edge for the latter case, a reflection of the complexities in numerics.
through which critical behavior (18) is introduced to the thermal properties of the system. In the original BEShydro Eq. (C8) is calculated and tested, but apparently Eq. (C9) is much more complicated to calculate numerically.
Here we validate the code by verifying that exactly the same results can be achieved using the two expressions for the Navier-Stokes term in Eqs. (C8,C9). A quick and easy test can be done using the setup from Sec. C 1, where are analytically available and easily used in BEShydro.
We have checked that both methods in Eqs. (C8,C9) give precisely the same results in Fig. 15. However, in the simulations of this work, the neos is used, and the two terms in Eq. (C10) are not analytical calculable. Thus we repeat the comparison to MUSIC as in Ref. [54] for the two gradients, and show the validation in Fig. 17, which show excellent agreement. From the figure, we also see that the baryon transport driven by D B ∇ µ n is very close to that driven by D B ∇ µ n+D T ∇ µ T , indicating that density gradients dominates over temperature gradients in the Navier-Stokes limit of the baryon diffusion current.