Dissipative hydrodynamics of relativistic shock waves in a Quark Gluon Plasma: comparing and benchmarking alternate numerical methods

This paper presents numerical cross-comparisons and benchmark results for two different kinetic numerical methods, capable of describing relativistic dissipative fluid dynamics in a wide range of kinematic regimes, typical of relevant physics applications, such as transport phenomena in quark-gluon plasmas. We refer to relativistic lattice Boltzmann versus Montecarlo Test-Particle methods. Lacking any realistic option for accurate validation vis-a-vis experimental data, we check the consistency of our results against established simulation packages available in the literature. We successfully cross-compare the results of the two aforementioned numerical approaches for momentum integrated quantities like the hydrostatic and dynamical pressure profiles, the collective flow and the heat flux. These results corroborate the confidence on the robustness and correctness of these computational methods and on the accurate calibration of their numerical parameters with respect to the physical transport coefficients. Our numerical results are made available as supplemental material, with the aim of establishing a reference benchmark for other numerical approaches.


I. INTRODUCTION
In its broadest and most general sense, hydrodynamics is the science of collective behaviour, i.e. the description of the dynamics emerging on top of the rules which govern the microscopic world. As a result, hydrodynamic behavior is expected to settle in the presence of a clearcut separation between the collective degrees of freedom and the microscopic ones from which they emerge. A safe margin in this respect is about two orders of magnitude, but the specific threshold at which such emergence takes place may vary from system to system and in some cases much less stringent thresholds may be present. Although usually associated with macroscopic motion, distinct hydrodynamic signatures can be found way deep into the microscopic world, down to the nanometric scale. In the last two decades relativistic hydrodynamics has also captured major interest from the apparently widely separate discipline of high-energy physics, the main driver being provided by the famous AdS/CFT theory, which sets an equivalence between d-dimensional field theory and (d+1) dimensional gravity [1]. This fascinating connection has opened an active field of modern research on so-called holographic fluids, i.e. strongly interacting quantum-relativistic fluids supporting the AdS/CFT duality. Among others, spectacular realizations of hydrodynamic holography have been reported for electron flows in graphene and quark-gluon plasmas [2][3][4]. Perhaps, most spectacular of all, recent experiments of the PHENIX collaboration [5], have reported evidence of the "smallest droplet ever", namely a droplet of quark-gluon plasma of the size of a few femtometers! These manifestations of hydrodynamic behaviour at truly short scales have spurred a major activity on the experimental, theoretical and, to a lesser extent, computational side. With specific reference to quark-gluon droplets, a major question pertains to the departures from strict hydrodynamic regimes which take place in small systems due to partial lack of "thermalization" of the initial strongly non-equilibrium configuration.
The quantitative description of such departures is beyond analytics, and raises major computational challenges, since a kinetic-theory treatment is not only computationally very demanding, but also theoretically questionable since strong-interactions imply supershort-lived quasi-particles, thus undermining the very premise of kinetic theory.
On the other hand, a macroscopic description of dissipative hydrodynamics poses further conceptual problems, since it has long been recognized that a naive relativistic extension of the Navier-Stokes equations is inconsistent with relativistic invariance, implying superluminal propagation, hence non-causal and unstable behavior. This can be corrected by resorting to fully-hyperbolic formulations of relativistic hydrodynamics. However, while various frameworks have been proposed, the definition of second-order relativistic viscous hydrodynamic equations is still debated, with a lot of ongoing research [6][7][8][9][10][11][12][13][14][15][16][17][18].
In this context, approaches based on a mesoscopic description help overcoming some of these problems. This approach is useful as, eventually, one may want to study "beyond-hydrodynamics" phenomenology; on the other hand, a consistent description of the hydrodynamic regime requires establishing a link between the mesoscale and the macroscopic parameters [10,12,14,[19][20][21][22][23][24][25][26][27][28][29][30][31]. More in general, there is a strong need for reliable numerical tools, as well as numerical benchmarks to compare accuracy, stability and performance of these solvers.
Typical benchmarks used in the validation of relativistic hydrodynamics codes are the Bjorken flow [32], Gubser flow [33], and the Riemann problem [34,35]. The latter is particularly useful to investigate the robustness of a numerical method due to the presence of strong discontinuities that give origin to the formation of shock waves. While the benchmark has an analytic solution only for two limiting cases, namely for the inviscid and the ballistic regimes, in the recent past several studies have analyzed the formation of relativistic shock waves in viscous QGP matter [36][37][38][39].
To the best of our knowledge, all previous works have been restricted solely to the study of the ultrarelativistic regime, for which the appropriate equation of state (EOS) writes as = 3P . In this work we consider instead a more general EOS, accounting for massive particles: where the relativistic parameter ζ = mc 2 /k B T , defined as the ratio between the rest mass energy mc 2 and the thermal energy k B T , physically characterizes the kinematic regime of the macroscopic fluid, with ζ → 0 in the ultra-relativistic regime and ζ → ∞ in the nonrelativistic one. We apply this benchmark and compare two different computational methods, and present numerical results exploring a wide range of parameters. We start by replicating results available in the literature in the ultra-relativistic limit and then move on to the study of fluids of massive particles, relevant to the QGP framework. The methods we consider both share a kinetic approach at the mesoscale level, but differ significantly in their numerical formulation, namely i) a Relativistic Lattice Boltzmann (RLBM) approach and ii) a Monte Carlo-enabled solution of the full kernel of the Boltzmann equation based on the Test-Particle-Method (RBM-TP). They are cross-validated for small values of the ratio between shear viscosity and the entropy density (η/s < 0.2 ), corresponding to a hydrodynamic regime where the Knudsen number is Kn 1. We show that the two solvers provide results in very good agreement, by analysing the profile of momentum integrated macroscopic quantities, as well as the non-equilibrium contributions to the moments of the particle distribution function [40,41].
We also show that, as expected, the RLBM approach fails when η/s is significantly large (ballistic limit) while successfully captures the physics features of flows at very low shear viscosity. The Monte Carlo approach, RBM-TP, is able to get solutions in agreement with RLBM even at very low viscosity, η/s 0.05, being able to naturally describe the evolution also for systems at large viscosity up to the ballistic limit η/s → ∞.
The numerical results presented in this work are made available as supplementary material with the aim of pro-moting further future comparisons [42].

II. MODEL EQUATIONS
The kinetic description of a relativistic gas is based on the particle distribution function f ((x α ), (p α )), depending on space-time coordinates (x α ) = (t, x) and momenta (p α ) = p 0 , p = p 2 + m 2 , p , with α = 0, 1, 2, 3. The space-time evolution of f ((x α ), (p α )) is governed by the relativistic Boltzmann equation which, in the absence of external forces, writes as C[f ] being the collisional operator. Macroscopic quantities are defined from the moments of the distribution functions. The first moment is the particle four-current: while the second moment defines the energy-momentum tensor: The balance equations of the particle four-current and of the energy-momentum tensor deliver the following conservation equations: These equations are purely formal until a specific form for N α and T αβ is specified. Following Landau and Lifshitz [43] one has with the energy density, P the hydrostatic pressure, q α the heat flux, π <αβ> the pressure deviator, the dynamic pressure, and ∆ αβ = U α U β −η αβ the (Minkowski-)orthogonal projector to the fluid velocity U α ; the latter, in the Landau frame, is defined as T αβ U β = U α .

III. NUMERICAL METHODS
In this section, we provide a brief account of the two numerical methods used in this work, the Relativistic Lattice Boltzmann Method and the Montecarlo Test-Particle Method.

A. Relativistic Lattice Boltzmann Method
The first computational method used in this comparison is a relativistic lattice Boltzmann method (RLBM). This method [44][45][46] is a computationally efficient approach to dissipative relativistic hydrodynamics. One key advantage over other relativistic hydrodynamic solvers is that, being based on a mesoscale approach, the emergence of viscous effects does not break relativistic invariance and causality, since space and time are treated on the same footing, i.e. both via first order derivatives (hyperbolic formulation).
This numerical method solves a minimal version of Eq. 2, where the microscopic momentum vector is discretized on a Cartesian grid, and with the collisional operator replaced by the Anderson-Witting single relaxation time approximation [47,48]: In the above, τ is the relaxation proper-time and f eq is the equilibrium distribution for which we consider the Maxwell-Jüttner statistics: here and in the following K i (ζ) is the modified Bessel function of the second kind of order i. The connection between the microscopic parameter τ and the macroscopic equations is given by the transport coefficients, for which we take into account the analytic expressions resulting from the first order Chapman-Enskog expansion [49]. Relevant for the present study is the shear viscosity η: with Ki 1 the Bickley-Naylor function

B. The Test Particle Method
The second computational approach considered in this paper also belongs to the realm of kinetic transport theory. We use a relativistic transport code developed to perform studies of the dynamics of heavy-ion collisions at both RHIC and LHC energies [21,[50][51][52][53][54][55][56][57].
Recently, the code has been extended to the solution of equation Eq. 2 for massive particles, which allows to simulate a fluid with an EOS that can be close to the recent lQCD calculations [57]. In this work, Eq. 1 has been used.
In this work we consider only 2 ↔ 2 collision processes, which give rise to a collisional operator C[f ] of the form: are the distribution functions of the outgoing (ingoing) particles, s = (p 1 + p 2 ) 2 and σ(s, Θ) is the differential cross section which is related to the total cross section by σ tot = dΩ dσ(s,Θ) dΩ . The numerical solution of the transport equation is obtained by using the test particle method, a popular option in many transport calculations [58][59][60]. In this method, the phase-space distribution function is sampled by mean of a large number of so-called test particles. In fact, it can be shown, that the phase space distribution given by a collection of point-like test particles is a solution of the Boltzmann equation, provided the positions and momenta of the test particles obey the relativistic Hamilton equations [61,62].
The numerical implementation of the collision integral is based on the so-called stochastic method [63] that has proven capable of describing efficiently also the ultrarelativistic limit, avoiding the issue of causality induced by a geometrical interpretation of the collision integral.
The transport code permits to study the dependence of physical observables on microscopical processes fixed by matrix elements or cross section. A novel approach is based on the idea of gauging the collision kernel C[f ] to a desired η/s ratio by an effective cross section σ tot [21,52,64] . Such an approach was inspired by the success of the hydrodynamical approach to describe experimental data [65][66][67] and permits to employ a Boltzmann transport equation in regimes where η/s (or equivalently the scattering relaxation time τ ∼ 1/σρ) is very low [54,68,69].
The expression given by the first order Chapman-Enskog expansion (Eq. 10), is used for the shear viscosity, with the relaxation time defined as τ = R −1 = 1 n σtrv rel , with σ tr the transport cross section, v rel is the relative velocity of the two incoming particles. In Ref. [21] it has been checked through the Green-Kubo correlator that employing Eq. 10 generates a fluid with the desired value of η/s.
For the massive case R is given by the following expression:

IV. NUMERICAL RESULTS
For our numerical analysis, we begin by considering the relativistic Riemann problem. This problem describes a tube filled with a gas which initially is in two different states on the two sides of a membrane placed at x = 0. As a result, the macroscopic quantities describing the fluid present a discontinuity at the membrane. Once the membrane is removed, the discontinuities decay producing shock/rarefaction waves, depending on the initial configuration chosen for the two different chambers. For simplicity we assume the flow to be homogeneous in the transverse directions.
Analytical solutions for this problem are available only for the massless case, and for extreme values of the ratio η/s, i.e. in the inviscid limit (η/s → 0) [34], and in the free-streaming limit (η/s → ∞) [70]. We use the same initial conditions as in [37], with a pressure jump defined as P 0 = 5.43 GeV fm −3 for x < 0 and P 1 = 0.33 GeV fm −3 for x > 0. Initial values for the temperature are respectively T 0 = 400 MeV and T 1 = 200 MeV.
We perform simulations at constant values of η/s. The calculation of η follows from the discussion in the previous section, while the entropy density is approximated using [49] we use the following expression for the equilibrium density with d G = 16 the degeneracy factor of gluons. Combining Eq. 10, 12 and 13 it is then possible to define the relaxation time required to keep the ratio η/s constant to a desired value k. As an example, the expression in the ultra-relativistic limit reads as: In RLBM simulations, the above equation is used to locally adjust the relaxation time, while in RBM-TP simulations it is used in combination with Eq. 11 to calculate the corresponding value of the transport cross-section.  For all values of the particle rest-mass m, we list the relaxation times corresponding to the initial conditions on the left (τ 0 ) and on the right hand side (τ 1 ) of the discontinuity.
As a warm-up exercise, we start by reproducing the results of previous studies in the ultra-relativistic regime. In particular, we compare against the Parton cascade Boltzmann Approach to Multi-Parton Scatterings (BAMPS) [71], which numerically solves the Boltzmann equation using a Monte-Carlo approach. In Fig. 1, we show that both the methods correctly reproduce the results provided by BAMPS at η/s = 0.1. RLBM also gives a good approximation to the analytical solution in the inviscid limit; in this case the test-particle methods cannot be applied since for η/s → 0 the cross section becomes unphysically large leading to numerical instabilities. On the other hand for large values of η/s, the hydrodynamic approach becomes questionable, as we transit towards a ballistic regime. Although at η/s = 0.5 RLBM still manages to capture the qualitative behaviour of the flow, in this regime the results provided by the RBM-TP are more reliable and in better agreement with BAMPS. It is however important to note that even at the very low η/s = 0.05, the RBM-TP is able to describe the dynamical evolution in excellent agreement with RLBM. Finally, in the free-streaming limit, the test-particle method reproduces correctly the analytic solution, while an unphysical "staircase" effect is observed in the profiles produced by RLBM; it has been shown that higher order schemes can cure this issue [72]. We remark that the RLBM simulations carry a systematic error due to the truncation of the higher order moments of the particle distribution (see [46] for details), whereas a statistical error is inherently associated with RBM-TP. However in this and in all other figures in this paper, error bars are not shown, since we average over a sufficiently large number of events in order to keep statistical errors well below 1% for the macroscopic observables.
In the following, we take into consideration fluids consisting of massive particles and restrict our analysis to values of η/s < 0.2, corresponding to a hydrodynamic regime for which we have observed a good agreement between the two methods. We recall that such a regime is also the one of interest for the study of the quark-gluon plasma (QGP) in ultra-relativistic collisions [73,74]. In Fig. 2, we show the time evolution of the pressure and velocity profiles for a specific case with m = 0.8 GeV and η/s = 0.1. We observe an excellent agreement in the dynamics predicted by the two solvers; in particular, we remark that the two cases compare well not only at times t much longer than the relaxation time τ , but also in the short-term regime, t ∼ τ .
In Fig. 4 we also show the space distribution of the relaxation time τ (x) in the two approaches for the case of η/s = 0.1. It is seen that, even if the initial values, dashed lines, are the same according to Eq. (10), the local evolution in the region of the shock wave exhibits different values for RLBM (solid lines) and RBM-TP (dots) for the massless case, a discrepancy that nearly disappears for the massive case m=2 GeV. This can be attributed to the different collision kernels: the Anderson-Witting in RLBM and the full Boltzmann in RBM-TP; however, as shown in all other results, such difference does not lead to any appreciable difference of macroscopic quantities, like the hydrostatic pressure or the collective flow. Some difference is observed instead in Fig. 6, which reports non-equilibrium quantities, such as the dynamical pressure or heat flux. Such differences also tend to vanish for the very massive case, as we are going to discuss at the end of this section. In Fig. 3, we show the results obtained considering again a rest mass m = 0.8 GeV but for different viscous regimes. The two methods are in good agreement, with only slight differences observed in the proximity of the shock-wave front. The curves for RBM-TP at η/s = 0.01 are not shown since, at such low viscosity, the methods becomes numerically unstable as the cross section and hence the computational time diverges as η/s → 0. However it has been shown in Ref. [54] that a linear extrapolation in 1/σ for σ → ∞ provides the correct pattern even for ideal hydrodynamics. In Fig. 5 we fix η/s = 0.1 and compare the results obtained for fluids of particles with rest mass of 0, 0.8, 2 and 4 GeV. Once again, both computational methods yield the same numerical results, which is remarkable given that no parameter fitting is performed in this analysis (apart from requirement to keep η/s constant).
Furthermore, it is interesting to analyse the nonequilibrium contributions to the four-flow tensor N α and to the energy-momentum tensor T αβ .
In the top panel in Fig. 6, we report the time and the spatial component of the heat flux, calculated via The time component of q α is expected to vanish for nonrelativistic flows. We can observe that as the rest mass is increased, the fluid moves at slower speeds and the correspondent peak in q 0 reduces accordingly. We observe significant discrepancies at the peak values of both the time and spatial components of the heat flux in the massless case. RBM-TP closely follows the results of BAMPS, while RLBM seems to be in better agreement with other lattice Boltzmann approach employing a collisional C[f ] kernel as the one in Eq. 8 for RLBM (see e.g. Fig. 6 in [72]). On the other hand, when considering massive particles the results seem to be in good agreement.
In the bottom left panel of Fig. 6, we show the T xx component of the energy-momentum tensor. The two methods provide results in good agreement, with the largest discrepancies observed once again in the ultrarelativistic limit.
Particularly relevant to the study of QGP is the analysis of effects related to the bulk viscosity, with potential implications for dark matter [75][76][77][78].
The bulk viscosity enters the non-equilibrium contribution to the energy-momentum tensor, and in particular is proportional to the dynamic pressure: We calculate the dynamic pressure from simulations by applying the projector ∆ αβ to Eq. 4, which yields: Measurements of are shown in the bottom right panel of Fig. 6. This quantity is exactly zero for the massless case, consistent with the fact that a ultra-relativistic gas has no bulk viscosity. Bulk effects are more prominent for the case m = 2 GeV, which corresponds to ζ ≈ 5, consistently with analytical calculations [49]. Differences in the dynamical pressure can have some relevance for the phenomenology of the QGP physics because they moderately affect the average momentum of the spectra and the build-up of azimuthal anisotropic collective flows [78,79].
In general, we can see that for the heat flux, q α , and for the dynamical pressure larger differences between the two methods, mainly in the region of the peak. This is not surprising, because such quantities focus on the details of the non-equilibrium dynamics and the simplified collisional kernel in RLBM may induce a slightly different non-equilibrium dynamics with respect to RBM-TP. However, this does not lead to any significant difference in the dynamic evolution of global quantities, like the pressure and the collective flow, as shown in all the previous section. It is also worth of note that the differences decrease with increasing mass, corresponding to a less rapid dynamics, and that the differences are larger in the same region where the τ (x), shown in Fig. 4, exhibits the largest departure between the two approaches. In any case, it should be observed that the differences between the two methods in are anyway smaller than 0.5%, with respect to the hydrostatic pressure.
We end this section by stressing again that our results, while conceptually based on a Boltzmann mesoscopic description, result from two fully different computational approaches; each computational method has its own parameters describing the mesoscale dynamics, but in both cases, these values can be derived from transport coefficients at the macroscopic scale. The two algorithms provide correct results in different ranges of the transport coefficients: RLBM works correctly in the range η/s → 0 to η/s 0.5, while RBM-TP works correctly for η/s ≥ 0.05 up to the ballistic limit of η/s → ∞. This suggests that the two methods could be used together to cover virtually the full range η/s parameter. However for the phenomenological study of the QGP dynamics in ultra-relativistic collisions both methods cover the relevant range usually explored, 0.08 ≤ η/s ≤ 0.2. The RBM-TP would be also the most suitable approach to study systems where the physical conditions are such that one evolves from the low viscosity to the ballistic regime, as can occur in electron fluids in graphene constrictions [80].
We finally mention that the RBM-TP provides a solution of the Boltzmann equation for the one-body distribution function f ((x α ), (p α )), allowing to evaluate not only the T µν components, but several physical quantities, like the transverse momentum spectrum, the azimuthal anisotropies v n as a function of the momentum, opening to a comparison with the wealth of experimental data. Furthermore in Ref. [54] it has been show that is possible to evaluate the viscous correction to the distribution function δf (p), an essential quantity for the solution of the Israel-Stewart viscous hydrodynamics. Of course, the access to such a wider class of observables comes at a price in computational time, that is currently more than two orders of magnitudes larger for RBM-TP w.r.t. RLBM.

V. CONCLUSIONS
Summarising, in this paper we have presented two numerical schemes for the simulation of relativistic hydrodynamics, namely the RLBM with the simplified Anderson-Witting kernel and RBM-TP that by test particles method solve the full Boltzmann collision kernel. While these two models build on different methodologies, it is quite remarkable that they yield to comparable results in the framework of relativistic hydrodynamics. In particular, a comparison of the two methods on a benchmark problem (the relativistic Riemann problem) shows matching results. The simulations have been performed at different values of the relativistic parameter ζ, and under different viscous regimes (i.e with different values of the parameter η/s).
In the ultra-relativistic case, we have correctly compared our results against the analytical solutions that are available in the inviscid and free-streaming regimes, as well as against data from a third numerical method, BAMPS, for all values of η/s between these two limits. It is worth noting that the test particle method is not suitable to simulate fluid flows in the limit η/s → 0, since in this region, numerical instabilities arise, due to unphysical values of the particles cross sections. However the comparison with RLBM shows that it is already very reliable even at very low viscosity such as η/s 0.05, smaller than the conjectured AdS/CFT lower bound. On the other hand, RLBM begins to struggle at values η/s > 0.5, where the hydrodynamic description starts to become questionable, although still capable of reproducing the general behaviour of the dynamics.
This suggests that the two methods could be applied in different zones of the parameter space, still featuring a window of cross-compatibility which permits their handshaking in prospective multiscale simulations of quarkgluon plasmas. The specific case of the phenomenology of QGP in ultra-relativistic collisions where η/s(T ) ranges in 0.08 < η/s < 0.2, is indeed within the range of crosscompatibility of the two methods.
After having established a firm connection to external results in the ultra-relativistic limit, we have tested the two numerical schemes on simulations at values of η/s, typical of quark gluon plasmas produced in colliders like RHIC and LHC and for several different values of the rest mass of the particles. The two methods are in excellent agreement when comparing the profiles of macroscopic quantities of interest, such as the hydrostatic pressure and the macroscopic velocity. Some limited difference emerges when comparing the non-equilibrium components of the particle four-current and the energymomentum tensor associated to heat-flux or dynamical pressure. For the future, would be interesting to explore whether the excellent agreement reported in this paper still holds for initial profiles closer to the initial stage of ultra-relativistic collisions, such as the Gubser [33] and the Bjorken profiles [32]. The results presented in this paper allow to benchmark the results of newly developed numerical tools against two different numerical schemes; this process of cross-validation would permit to test the accuracy of such schemes, thus paving the way to reliable and accurate simulations of real physical systems, such as quark-gluon plasmas dynamics in current and future high-energy experiments.
ACKNOWLEDGMENT DS has been supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 765048. SS acknowledges funding from the European Research Council under the European Union's Horizon 2020 framework programme (No. P/2014-2020)/ERC Grant Agreement No. 739964 (COPMAT). Numerical work has been performed on the COKA computing cluster at Università di Ferrara, and the computing cluster GR4-LNS at the INFN-LNS and the QGPDyn cluster at the University of Catania.