Inclusive and effective bulk viscosities in the hadron gas

We estimate the temperature dependence of the bulk viscosity in a relativistic hadron gas. Employing the Green-Kubo formalism in the SMASH (Simulating Many Accelerated Strongly-interacting Hadrons) transport approach, we study different hadronic systems in increasing order of complexity. We analyze the (in)validity of the single exponential relaxation ansatz for the bulk-channel correlation function and the strong influence of the resonances and their lifetimes. We discuss the difference between the inclusive bulk viscosity of an equilibrated, long-lived system, and the effective bulk viscosity of a short-lived mixture like the hadronic phase of relativistic heavy-ion collisions, where the processes whose inverse relaxation rate are larger than the fireball duration are excluded from the analysis. This clarifies the differences between previous approaches which computed the bulk viscosity including/excluding the very slow processes in the hadron gas. We compare our final results with previous hadron gas calculations and confirm a decreasing trend of the inclusive bulk viscosity over entropy density as temperature increases, whereas the effective bulk viscosity to entropy ratio, while being lower than the inclusive one, shows no strong dependence to temperature.


I. INTRODUCTION
Transport coefficients give insights about the microscopic dynamics of interacting matter close to equilibrium. The shear viscosity over entropy density η/s is the most extensively studied transport coefficient in relativistic heavy-ion collisions (RHICs); since the first viscous hydrodynamic calculations became available in 2008 [1,2], the extraction of η/s and its temperature dependence has been increasingly refined over the last decade. However, the situation is a bit different for the case of the bulk viscosity. Since AdS/CFT calculations (as models for QCD dynamics in the strong coupling) imply that it is very small for nearly-conformal systems [3], the bulk viscosity ζ (and its corresponding dimensionless ratio ζ/s)-which can be thought of as the resistance to uniform expansion/compression of a fluid-has not been subject to the same extended treatment as η in the context of RHICs [4]. It should be pointed out that although ζ is identically zero in conformal fluids [5], and that QCD approaches conformality in the limit of high energies/temperatures [6,7], there is no evidence that the nuclear matter which is produced in accelerators (even at the highest LHC energies) is a conformal fluid. Moreover, the system becomes less and less scale invariant as the system cools down with time [8].
Although not exhaustive, some studies on the effect of bulk viscosity on some observables such as elliptic flow [9, 10] and particle spectra [11] were made. More recently, the bulk viscosity has started attracting more attention since it was pointed out by phenomenological studies in hybrid models that the inclusion of bulk viscosity as described by [12] was important in some cases to properly reproduce simultaneously the radial and azimuthal flow anisotropies [13,14]. Most notably, the first quantitative extractions of shear and bulk viscosities employing Bayesian techniques have recently appeared [15,16]. In these works, the functional form of the temperature dependence of the transport coefficients influences the prior and therefore an external input for these is very important. In particular, the bulk viscosity is expected to have a peak around the transition from hadronic matter to the quark-gluon plasma [8]. Close to the phase transition at vanishing baryochemical potential, lattice QCD calculations indicate a large enhancement of the bulk viscosity [8,17]. Above the crossover temperature, a fast drop-off is also suggested by quasi-particle models [18,19].
On the purely hadronic side, theoretical calculations of the bulk viscosity are notoriously more complicated than those of the shear viscosity, and as such are scarcer. However, using different models and computational techniques, the temperature dependence of the bulk viscosity of a hadron gas was presented e.g. in Refs. [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. Results from the various calculations differ from one another by an order of magnitude or more, as we will see when comparing our own results with some of these calculations.
Among these calculations we will pay special attention to those restricted to very low temperatures where pions dominate the hadronic mixture. In this regime the interactions of pions can be described by chiral perturbation theory (and its unitarized version to describe the resonant energy domain). One of these calculations [24] applied a diagrammatic Green-Kubo method and predicted a double bump structure for ζ at low T . The first of these bumps was explained from the explicit conformal breaking due to the pion mass, while the second was related to the conformal anomaly appearing at temperatures close to the crossover. A calculation with similar interactions but using the Boltzmann-Uehling-Uhlenbeck equation [29] further commented that the addition of a pion pseudochemical potential was also necessary for a consistent treatment of pion elastic collisions. However, in Ref. [28] the focus was on the much slower 2 ↔ 4 pion inelastic processes and obtained a very different value of the bulk viscosity (diverging at T = 0). These rather different calculations illustrate the effect of including or excluding particle number-changing processes: whereas Ref. [28] uses the idea that the slowest processes (inelastic collisions) should dominate the value of ζ, Refs. [24,29] argue that such processes are so slow that they cannot be effective at all in RHICs. In this paper we will clarify the conceptual difference between the two points of view-distinguishing between "inclusive" and "effective" bulk viscosities-by addressing this coefficient using a microscopic simulation code, SMASH (Simulating Many Accelerated Strongly-interacting Hadrons).
In the following we will present various results for the bulk viscosity in simple hadronic systems of various chemical compositions, amongst which hadronic predictions for hydrodynamical calculations of RHICs. Some more technical considerations that have to be taken into account in order to obtain them will also be discussed.
In Sec. II we introduce the methodology to extract the bulk viscosity via a Green-Kubo relation and introduce the SMASH transport approach. In Sec. III we apply the model to a simple relativistic gas interacting with constant cross section, where comparison with the corresponding Chapman-Enskog solution will calibrate our model in terms of systematic uncertainties. In Sec. IV we show that adding resonances to the system requires a revisiting of the assumption made for the form of the correlation function. We show how the simple exponential decay ansatz breaks down, and further analyze the effect of the resonance lifetimes. In Sec. V we apply the method to the full hadron gas for several temperatures and box sizes. We introduce definitions for the inclusive and effective bulk viscosities and present final results for both ζ/s and ζ eff /s, comparing with previous calculations. Finally, in Sec. VI we summarize our work.

II. METHODOLOGY
A. Green-Kubo formalism In this work we apply the Green-Kubo formalism [36][37][38] to obtain the bulk viscosity coefficient of different systems. While different versions of the Green-Kubo formula exist in the literature depending on the system and thermodynamical ensemble used, the most general form reads [39][40][41][42][43] where V is the volume of the system, T is the temperature, and ∆Π(t) ≡ Π(t) − Π is a fluctuation around the thermodynamical equilibrium average. The variable Π is defined as where P (t) = 1 3 T i i (t) 1 is the (instantaneous) pressure, (t) = T 00 (t) the (instantaneous) energy density and n(t) = j 0 (t) the (instantaneous) particle density. All components of the energy-momentum tensor T µν and the particle 4-current j ν are understood to be averaged over V , In Π(t) appear two thermodynamical quantities: the speed of sound at constant number density and the compressibility at constant energy density. These quantities naturally appear in the source function (left-hand side) of the Boltzmann equation when considering the bulk viscosity of a gas with a conserved (net) particle number [20,[44][45][46].
For later reference the adiabatic speed of sound at constant entropy per particle S = s/n is related to these two quantities as [44,46] where w = + P is the enthalpy density. Expressions for all these quantities as functions of temperature are given in App. A. For convenience let us define the autocorrelation function which will be extracted from our numerical SMASH simulations and integrated over time as in Eq. (1). For other transport coefficients such as the shear viscosity or electric conductivity in dilute systems [47][48][49][50][51][52] it is generally assumed that the correlation function takes the form of a decaying exponential. This ansatz can be motivated by the relaxation-time approximation of the Boltzmann equation [53] or by the causal hydrodynamic equations [54]. This assumption should always be checked a posteriori within the precision of the data acquired. For the bulk viscosity, it is thus assumed that where τ ζ is the bulk relaxation time of the system. From Eq. (1), it follows that In some previous works the relaxation time τ ζ has been estimated to be related to the mean free time of the particles, i.e. the average time between collisions used for other transport coefficients as well. However this introduces a new source of uncertainty, as different transport coefficients are sensitive to different transport mechanisms. For example, while the mean free path is inversely proportional to the total cross section, the shear viscosity is sensitive to the "transport cross section", which could be a factor of 2 smaller than the total cross section in a p−wave scattering [50]. More importantly, using the mean free path misses the dependence of τ ζ on the resonance lifetimes, which was noted to be very important in the shear viscosity case [50]. As we will see, this will also prove to be particularly significant for the bulk viscosity.
It is helpful to realize that the value of C ζ (0) is an equilibrium quantity. From its definition, where f eq 1 = f eq (p 1 ) is the (spatially-averaged) distribution function in equilibrium e.g. the Maxwell-Boltzmann function f eq 1 = g exp[−(E 1 − µ)/T ] (g is the internal degeneracy of the particle), and E 1 = p 2 1 + m 2 . To compute C ζ (0) we need to know the equal-time correlation function of the spatially-averaged fluctuation of the distribution function, for which we can directly apply the result of [55] for the 2-point correlation function, Combining these, we obtain Incidentally, this formula exactly coincides with the quantity T ζ/τ ζ R derived in [46] for a (Bose) gas with binary interactions. The result in [46] uses the relaxation time approximation, where one identifies τ ζ R τ ζ . Using the expressions given in App. A for the different thermodynamics quantities appearing in (11), one can compute the explicit temperature dependence of C ζ (0).

B. Hadron gas modeling: SMASH
In this work we use the SMASH transport approach [56,57] to simulate infinite hadronic matter in a box with periodic boundary conditions. In SMASH, all well-established hadrons of the PDG [58] are included, with their interactions modeled by resonance excitation and decay, elastic as well as inelastic 2 ↔ 2 processes.
At this point it is important to mention that the V used in our analysis is the entire simulation box volume, instead of a subvolume of the whole system. By doing so, we get that the total energy and total particle number are conserved (at least in simple systems), bringing the system in a sort of microcanonical ensemble over V . Therefore ∆n(t) = ∆ (t) = 0 , and the correlation function reduces to The instantaneous pressure P (t) is extracted from the energy-momentum tensor T µν (t) of the equilibrated system, following the methodology described in [50,51]. Such simulations provide the complete phase-space information of all particles in the system, which are in this case discrete, and given at specific time steps. For this situation, we can define the components of the energy-momentum tensor as where N is the total number of particles in V , p µ i is a component of the momentum 4-vector associated with particle i.
The averaging contained in the correlation function (12) also has to be defined for the discrete times t ≡ u∆t at which the information is available, where K is the total number of considered time steps, u is a positive integer with u < K and ∆t is the time interval between each time step. It is numerically challenging to take the limit of K → ∞ in Eq. (14) and thus the relative error of any numerical computation of the correlation function necessarily increases rather quickly with time and eventually reaches a state of pure noise, as one can see for example on Fig. 2.

III. SIMPLE GAS WITH ELASTIC INTERACTION
A single-component relativistic gas interacting through elastic collisions provides the first example to test our method. In the case of a gas with constant, isotropic cross section (hard-sphere gas) the bulk viscosity is zero in the nonrelativistic and the ultrarelativistic limits [5,8]. However in an intermediate domain of temperatures the bulk viscosity is small, but nonzero. Without loss of generality we will assign a mass to the particles m = 138 MeV, and internal degeneracy of g = 3 (resembling a pion gas but interacting with a constant cross section of σ = 20 mb).
Such a hard-sphere gas has been studied before in the context of the bulk viscosity. Its value has been extracted analytically e.g. in [59] by linearizing the collision term of the Boltzmann equation using the Chapman-Enskog approximation to first order (see also [32]). More generally, by modifying the numerical codes used in [29,46] we can extend the Chapman-Enskog expansion for this system to higher orders to check convergence. This will help us to calibrate the Green-Kubo calculation in this simple case. To start with, it is instructive to look at a sample of the measured fluctuations of the pressure in such a system and to compare it to the off-diagonal energy-momentum fluctuations T xy (t) associated with shear viscosity (see [50,60]), which is done on Fig. 1. While the fluctuating nature of the two signals appears relatively similar at first glance, there are significant differences between them. First, the amplitude of the signal for T xy (t) is roughly 25 times larger than the case of the pressure (notice the different OY axis scales). Second is the fact that pressure does not oscillate around zero, and thus an average pressure needs to be subtracted in the correlation function. This is not as trivial as one could think, as the average pressure also introduces a statistical error which can be non-negligible. While the calculation of the correlation function is done over 4000 time steps spanning 200 fm as in the case of the shear viscosity (see [50]), we find that in order to get results in which the statistical error does not completely wash out the signal, the averaging of the pressure requires much larger data sets. We determined that for the studied cases, an averaging going over 100 000 time steps spanning 5000 fm was sufficient, in the middle of which we perform the previously mentioned calculation of the correlation function.
Note that, in principle, the thermodynamic pressure can be calculated analytically for such a gas assuming Boltzmann thermodynamics (e.g. via the J n,k functions defined in App. A). However, in more complex systems, although the SMASH equilibrium is very close to the grand canonical one, it can deviate from it slightly. For these bulk viscosity calculations, even minimal deviations in the average pressure of the order of a fraction of a percentage point can make a significant difference in the final signal, and thus such an analytical calculation would be highly non-trivial to perform to the required precision. Therefore, to keep our methodology consistent with the following sections, we always use the numerical extraction of the average pressure, instead of the Boltzmann expression.  Figure 2 shows a collection of correlation functions. The left panel illustrates how a rising temperature leads to a steeper C ζ (t) (which translates to a shorter relaxation time) as well as the expected increase of the statistical error as time increases. What is quite unique to the case of the bulk viscosity is that the initial value C ζ (0) can have a relatively large statistical error, up to 20% in this case, whereas in previous works the same error on the shear viscosity [50] or electrical conductivity [51] correlation function initial value was never larger than ∼ 6%, which would barely be visible.
The right panel of Fig. 2 additionally shows that the size of the box used for the calculation scales as 1/V [cf. Eq. (11)], so that the increase of factor of 3 in the size of the box is reflected by a decreasing of a factor of 27 in C ζ (t), with the slope being the same in both curves. Reducing the size of the box also reduces the relative statistical error, as the size of the fluctuations with respect to the average pressure diminishes as volume increases (10).
To calculate ζ, one then has to strike a balance between having a system which is large enough for thermodynamic quantities to be calculated, but small enough that the signal does not get washed out by the statistical error. This volume might differ for each value of the temperature within the same system, as can be seen in Table I, where we provide the specific volumes used for each temperature. Eq. (11) at the corresponding temperature. Larger/smaller boxes refer to the values presented in Table I; notice how the bigger volumes provide larger error bars, which is consistent with the larger error on the correlation function for these volumes. We observe that a very good agreement is obtained between the two, providing a nontrivial check on the method.
We proceed to fit the correlation function to the exponential decay form (6). Notice that the relatively large uncertainty in C ζ (t) makes it difficult to systematically decide where to stop the fit; we will simply stop it at t = 5 fm for this simple gas case.  Table I. The theoretical calculation comes from the Chapman-Enskog estimates in Refs. [59] (1st order) and [46] (first to third order). Right panel: Bulk relaxation time for the same system as a function of the temperature, for the same box sizes.
The bulk viscosity of this system calculated using the Green-Kubo formalism is compared to the semi-analytic Chapman-Enskog expansion in the left panel of Fig. 4 as a function of the temperature. The first-order Chapman-Enskog result is taken from [59] and numerically re-calculated with the method of [46], which also allows to go to third-order Chapman-Enskog where convergence is achieved. The agreement is rather good for temperatures between 100 and 175 MeV, even for smaller system sizes. At low temperatures, the agreement starts to break down, and, although not shown here, crumbles completely at even lower temperatures. At those low temperatures we observe that the correlation function is still exponential, but the uncertainties are large: the number of pions at these temperatures is so small, that statistics are very poor, making results at lower temperatures than shown unreliable. In parallel, using large volumes to increase statistics also washes out the signal; this can be seen in the black dots, which not only underestimate the analytic ζ but also see their error bars increase significantly. We are thus unfortunately not able to observe the nonrelativistic limit in which the bulk viscosity turns to zero at T → 0. In the right panel of Fig. 4 we show the relaxation time τ ζ . The uncertainty of the large volumes is again much larger than the smaller volumes, although the central values are compatible with those of the smaller volumes. The two panels of the figure show a good agreement between different calculations, validating the method for more complex systems.

IV. EFFECT OF RESONANCES
It is known that the presence of internal dynamical degrees of freedom (rotational, vibrational) as well as inelastic collisions-allowing a redistribution of internal energy in a more efficient way-contribute critically to the bulk viscosity [61,62]. The latter might happen via strongly number-changing processes like the 2π ↔ 4π considered in [28] or the NN → 5π annihilation, but also due to the presence of continuous resonance decay and recombination. These processes made a notable difference already for the shear viscosity [60], and their role is expected to be even more relevant due to the nature of the bulk viscosity coefficient.
We start this discussion by presenting the bulk correlation function for the full hadron gas with resonances. While we relegate its full analysis to Sec. V, it will first serve us as a motivation for the consideration of a more general ansatz for C ζ (t) in the presence of several hadron species and resonances.
A. Breakdown of the single exponential ansatz A solid baseline has been established for the calculation of the bulk viscosity at temperatures between 100 MeV and 175 MeV after using a simple pion gas with constant cross section. We directly proceed to the case of the full hadron gas as described by SMASH v1.6 [57]. As mentioned earlier, this gas includes not only elastic but also inelastic processes, be they binary inelastic 2 ↔ 2 interactions or, most commonly, indirect resonant 2 ↔ 1 ↔ 2 reactions where two particles will form a resonance of mass m, width Γ(m) and a sampled lifetime averaging at τ life = 1/Γ(m), after which it will decay into two new daughter particles which can or not be of the same species as the original ones (it is also possible for resonances to scatter and form larger resonances with other particles during their lifetime; see [56] for details). Note that in order to calculate the average pressure of this system to an appropriate degree of precision, we require simulations to provide at least 5000 fm of equilibrium data; this is extremely costly in terms of computational power, and as such limits the breadth of the exploration of the parameter space. While in this system the particle number is not formally conserved by the inelastic processes, we make the approximation that the contribution of the third term in Eq. (2) is small with respect to the pressure fluctuations. We checked that particle number fluctuations are of the same order of magnitude as pressure fluctuations, but the former are multiplied by a small (∂P/∂n) (see Fig. 12), largely reducing their contribution.
Let us consider the normalized (i.e. divided by their value at t = 0) correlation functions presented on Fig. 5 for different temperatures. As is readily visible, these offer a considerably different picture as what we observed in the previous case in Fig. 2. First, the statistical errors are here much less significant than they previously were. This is expected, as the introduction of resonances (and thus of mass-changing processes which dissipates the otherwise purely kinetic energy) leads to a massive increase in the magnitude of the fluctuations with respect to the average pressure, and as such, it is expected that the error on the pressure plays a smaller role in this case.
More importantly, the correlation functions at all temperatures display a somewhat peculiar shape. In the first 2-3 fm a period of rapid exponential decorrelation is followed later on by a less abrupt decay over relatively long times before the relative error finally increases to a point where the signal is dominated by noise. It is evident that the correlation functions are not describable by a single exponential function, and one needs to abandon the simple ansatz in Eq. (6). To physically understand these important modifications of the shape of the correlation function, we look at a toy system with a minimal content of particles. Let us consider a box with pions (with their physical mass and isospin degeneracy) interacting through a single resonance, the ρ meson. We switch off all other possible resonances and set to zero any contact cross section. This is a relatively simple system, in which we scale the lifetime of the decay ρ → π + π by a multiplicative factor. FIG. 6. Bulk correlation function at T = 150 MeV for a pion gas using the ρ resonance for the cross-section for the cases where the ρ has zero, a fifth, half or its full lifetime. Figure 6 shows the correlation function at T = 150 MeV of the π − ρ mixture in which the lifetime (and thus its relative abundance 2 ) is varied and eventually taken to zero. In this precise limit we recover the previously discussed case of 2 ↔ 2 elastic scattering with no intermediate resonance. First notice that the correlation decays exponentially, as we previously assumed. Second, it is also evident that even a very small τ life profoundly modifies the underlying physics. Such a lifetime allows for a continuous formation and decay of a resonance allowing the imbalance in the longitudinal (bulk) channel to relax in a more effective way than a pure local collision. Such a mechanism produces a large increase of the fluctuations (seen in C ζ (0)), and a reduction of the relaxation time. By decreasing the resonance lifetime, we increase the number of decays and recombinations per unit time (π + π → ρ → π + π). Therefore the relaxation time of the bulk viscosity is shortened, as can be seen in the figure. However, in the zero lifetime limit the collisions become effectively elastic, and we reach a limit in which the momentum relaxation is ineffective, with a very large relaxation time. This example shows the large dependence of τ ζ on the resonance lifetime. However, no deviation from the exponential form can be inferred so far.

C. Several resonances: effect on correlation function
Our previous analysis concerning the relaxation time dependence on the resonance lifetime was still possible on the basis of the single exponential decay of C ζ (t). For such a system with a single channel (one resonance) the correlation function does not develop a non-exponential behavior similar to what we see in Fig. 5. As the next step, it is possible to speculate that the presence of several interactions and decay modes, with a variety of relaxation times, determines the more complicated form of C ζ (t).
To validate this hypothesis we study the effect of two independent resonant channels in the box. We start with the π − ρ system of the previous section (using the physical ρ lifetime, not modified anymore). In addition, we introduce a parallel particle/resonance system in the simulation. We add a non-physical particle species B with the same mass as the pion (m = 138 MeV) interacting through a single resonance B * with the same pole mass as the ρ meson (m = 776 MeV); however, we use a much smaller decay width for the B * (see Table II). In summary, we insert a copy of the π − ρ system but with a longer-lived resonance 3 . Finally, to simplify the analysis, note that the π − ρ and the B − B * are not coupled to each other. The reason for such a particular system is the following. Under a fluctuation in the bulk channel, the π−ρ subsystem will have a relaxation time of the order of τ life ∼ 1/Γ 1 fm (similar to the result in the previous section). The new B − B * system has a significantly lower cross-section, and a lifetime which is an order of magnitude larger; this new subsystem is thus expected to relax ∼ 10 times slower than the π − ρ one, and it is expected that this separation of time scales will be visible in the correlation function of the mixture.
We plot the bulk correlation function of the different systems in the left panel of Fig. 7. As expected, the π − ρ subsystem (in blue) has a smaller relaxation time than the B − B * system (flatter red curve). The smaller C ζ (0) of the B − B * is due to the more suppressed resonant contribution, as the broader ρ resonance weights more in the thermodynamic average of C ζ (0). For both subsystems the correlation function is a single exponential, as expected.
Looking at the correlation function of the mixture of π −ρ−B −B * , we observe a non-exponential shape comparable to the ones for the full hadron gas in Fig. 5. Even more interestingly, adding up the individual exponential contributions from the π − ρ and B − B * systems results very precisely in the same correlation function for the full system, with later times being dominated by the B − B * process, the slowest one. We proceed to fit the resulting correlation function of the mixture at T = 150 MeV to a double exponential form, The tail of the correlation function is first fitted to extract C ζ,B (0) and τ ζ,B , and then subtracted from the total correlation function. After checking that the remaining function is indeed exponential, it is fitted to obtain C ζ,π (0) and τ ζ,π . The final fit is shown in the right panel of Fig. 7 in dashed line (the correlation function itself is hidden by the fit, but its error band is still visible). We obtain τ ζ,π = 0.91 fm and τ ζ,B = 9.65 fm. These values turn out to be of the same order of magnitude as the respective lifetimes 1/Γ ρ 1.32 fm and 1/Γ B * 9.85 fm. Thus, both microscopic processes of resonance formation/decay do affect the bulk viscosity of the mixture, each of them with its own characteristic relaxation time.
It is therefore natural to expect that the full hadron gas, being a massively more complex system, will be described by a collection of individual exponentials. However, contrarily to the case we just discussed, since many of the subsystems of the full hadron gas are actually coupled to each other, it would be complicated to associate these individual exponentials to a specific individual particle-resonance pair. Each one will correspond to each of the many interlinked subsystems (containing elastic and/or inelastic processes) present in the gas, with later times being dominated by the slowest such subsystem (i.e. the one containing the slowest set of processes).
In a more general way, one should then replace the single exponential ansatz by a linear combination of many such exponentials, where the kernel function of relaxation times ρ(τ ) is normalized to ∞ 0 dτ ρ(τ ) = C ζ (0) and can be found, in principle, via inverse Laplace transform of the correlation function [63]. Notice that the range of the possible relaxation times runs from 0 to +∞, accommodating fast as well as slow processes.
However, in the remaining part of this work we do not need to use the full integral version of Eq. (16), as we will see that the kernel function ρ(τ ) can be taken as a linear combination of a few Dirac deltas, one for each relaxation time taking place in the system. Notice that for N = 1 one recovers Eq. (6).

V. FULL HADRON GAS
We focus again on the correlation functions of Fig. 5 for the full hadron gas, and use the multi-exponential form (16) and (17) to fit them. By inspection, we observe that N = 3 components (that is, three Dirac deltas) are sufficient to achieve a good fit of the correlation functions. The corresponding relaxation times should be considered as the dominant modes contained in the kernel ρ(τ ), which are related to physical processes in the hadron system. Of course, many other relaxation times do exist in the mixture (in principle, as many as independent microscopic processes), but they carry such a small amplitude that are not reflected in the correlation function. We perform a global fit using the ROOT library, which takes into account the error band of C ζ (t) and provides statistical uncertainties of the fitting parameters. A much more detailed discussion on the multi-exponential fitting can be found in App. B.
In Fig. 8 we plot the final result of the bulk viscosity (left panel) and the bulk viscosity over entropy density (right panel) for the full hadron gas as functions of the temperature. We provide the result for two different sizes of the box, denoted as "larger" and "smaller box", the precise lengths of which are given in Table III Table III.
Along with this "inclusive" ζ and ζ/s (where all modes present in the fit are included in the calculation) we have included results of an "effective" bulk viscosity coefficient. The latter is calculated by taking the long-lived modes out of the analysis for phenomenological reasons: In an infinite-lived system all modes contribute to the correlation function at some point, as the total relaxation of a fluctuation does not happen entirely until all modes in the system have equilibrated. In particular, the slowest mode is typically the one that dominates the bulk viscosity, as it is the one describing the long tail of the correlation function; however such slow processes are not effective in a short-lived system, if their inverse rate is much larger than the lifetime of the system. If in RHICs the hadronic phase lasts approximately 10-30 fm/c, then a relaxation mode with τ ζ = 10 2 − 10 3 fm cannot play any role. The part of the system corresponding to that mode remains out of equilibrium for the whole time, and does not contribute to the transport coefficient calculation. We define the effective bulk viscosity ζ eff to be the transport coefficient where such modes have been excluded.
More formally, the effective bulk viscosity can be defined as where the effective correlation function now depends on a cutoff τ * , or the order of the lifetime of the system, above which the modes are suppressed. Using e.g. a hard cutoff to remove these modes, Note that to obtain the effective bulk viscosity one still integrates the correlation function up to ∞, but the kernel ρ(τ ) is restricted. This definition still assumes the validity of the exponential ansatz for every mode. Why should this effective bulk viscosity be of any relevance? Suppose that one tries to describe the evolution of the system by a relativistic hydrodynamic code for heavy-ion collisions with the bulk viscosity as an input parameter to be fixed 4 . We argue that the extremely long-lived processes will hardly happen during the real evolution of the system, so they cannot be part of the eventually-inferred viscosity. The effective transport coefficient defined here should be associated to the one obtained from matching experimental observables using hydrodynamic codes; in contrast, the inclusive bulk viscosity should rather be compared with a theoretical calculation, e.g. solving the Boltzmann equation in the thermodynamic limit. Due to the suppression of the dominant mode (or modes), it is clear that ζ eff should always be smaller than ζ.
To completely understand this distinction, let us finally present another example of such an effective viscosity. Mannarelli et al. [64] calculated the shear viscosity due to phonons in optically trapped cold Fermi atoms. At low temperatures, the mean free path of phonons increases and exceeds the physical boundaries of the superfluid region. The shear viscosity is proportional to the mean free path, so at low temperatures it is possible to define an effective shear viscosity where the mean free path is replaced by a distance of the order of the atomic cloud. In our particular case, the bulk viscosity is proportional to a linear combination of relaxation times, and the effective bulk viscosity imposes a limiting time of the order of the system's own duration. In our study, motivated by RHIC physics, we calculate this effective bulk viscosity by removing the slowest mode of the three (last column in Table IV), whose relaxation time is typically much larger than the hadronic lifetime in a RHIC. Notice that for the highest temperatures T = 142 MeV and T = 172 MeV, τ ζ,3 is actually of the order of the lifetime of the fireball, and one could argue that this mode can still play some role in heavy-ions. Therefore, one should strictly interpret the effective bulk viscosity as a lower bound in these cases. Also note that this implies that systems with different lifetimes could have a different effective bulk viscosity, such as for example in the experiments at the very different beam energies of the Relativistic Heavy Ion Collider and the Large Hadron Collider.
The final results for ζ and ζ/s in Fig. 8 behave similarly for both box volumes. ζ/s decreases systematically with temperature to values around ζ/s 1, perhaps reaching a plateau around T = 172 MeV. Only the result for T = 86 MeV is quite different in the two volumes, and might correspond to a poor quality in one of the volumes used, similar to the discrepancy in the simple pion gas in Sec. III.
The effective bulk viscosity is always smaller in magnitude, as expected. The ζ eff at T = 172 MeV is somewhat different between the two volumes, due to the different value of the τ ζ,2 for that temperature. It is not evident to us which one of the two, if any, is of lesser quality. ζ eff /s is a rather flat or slightly increasing function of the temperature for the considered range.
The final value of our coefficients is obtained by averaging the two box sizes and combining their uncertainties. Our average value for ζ(ζ/s) and ζ eff (ζ eff /s) is shown in the left (right) panel of Fig. 9.

A. Discussion and comparison
In this section we attempt to contextualize the present calculation by testing it against previous calculations of the bulk viscosity. Before doing so, let us briefly comment on the relation between the bulk viscosity and the adiabatic speed of sound v S defined in Eq. (4).
In a massless, weakly-coupled gas, previous calculations using the Boltzmann equation and kinetic theory have shown that the relation between shear and bulk viscosity should be proportional to the squared non-conformality parameter [43,61,65,66] [67], and the lattice/QCD data is extracted from Ref. [6].
We can try to estimate the adiabatic speed of sound from the measurements of the energy density and pressure in SMASH for the full hadron gas studied before. This is illustrated in the left panel of Fig. 10, where we plot the values of these two quantities for each of the four temperatures. This shows the dependence of P versus needed to obtain the speed of sound. Before extracting v 2 S , we verify whether the entropy density s, number density n or entropy per particle s/n is held constant in these measurements, as we have not imposed any of those conditions explicitly. This is detailed in Table V, where we provide the values of s, n and s/n for each temperature. None of the three quantities remains absolutely constant but it is clear that one can rule out an isentropic (constant s) and isochoric (constant n) dependence. On the other hand, the entropy per particle does not vary much. Therefore it is fair to assume that the speed of sound, obtained from the relation between P and in our plot, will be approximately adiabatic (constant s/n), or, at least, a close proxy for it. We parametrize the dependence of the pressure to the energy density with a power law, and find that P ( ) = 0.153 0.914 , where both quantities are measured in GeV/fm 3 . The fit is shown as a solid line in the left panel of Fig. 10. In the right panel of the same figure we show the resulting v 2 S from this relation in blue dots, which is a decreasing function within this range of temperatures. Our values compare well with the result of Ref. [67] for a hadron resonance gas including resonances up to a mass of 2.5 GeV (similar to ours), and it is also comparable with the lattice QCD calculation of Ref. [6], the deviation at high temperatures being due to the absence of a deconfined phase in our model. Finally we move to a comparison of available calculations for ζ/s presented in Fig. 11, and we shortly discuss every other result with our own. We have shown both ζ/s and ζ eff /s computed from SMASH in red and green symbols, respectively.
• Noronha-Hostler et al. [23] use a hadron resonance gas model which assumes a comparable set of hadronic states as the ones used by SMASH. However, this model assumes a noninteracting tower of states, and the hadron resonance gas is supplemented with an exponentially increasing density of Hagedorn states. The bulk viscosity is calculated using the small-frequency spectral ansatz presented in [8], which matches the Euclidean version of the correlator of the trace of the energy-momentum tensor. Their result is comparable to our ζ eff /s, as that calculation lacks of the very slow dynamical process affecting our viscosity. An increase of ζ/s close to T c is only obtained by the inclusion of the Hagedorn states, and such an increase is not captured by any other model, except perhaps, by the SMASH effective bulk viscosity.
• The Rougemont et al. calculation is performed in a holographic setup, and as such is difficult to compare to our own results. The bulk viscosity is small in this calculation and rather flat, which is compatible with our value of ζ eff /s.
• The calculation of Dobado et al. [29] and Lu et al. [28] are both computed for a pure pion gas using chiral perturbation theory at low temperatures. However, their different approaches illustrate the conceptual difference between ζ and ζ eff . While [28] considers the slowest number changing process affecting the bulk viscosity (2π ↔ 4π) and neglects any elastic collisions, the calculation in [29] does not consider this process and uses the 2π ↔ 2π process only with a pion pseudochemical potential. In the first calculation the extremely slow inelastic process (suppressed by the derivative coupling at low energies) describes ζ. In the second calculation these processes are absent during the hadronic stage of RHICs and only elastic collisions are able to build a ζ at the expense of the change in chemical potential. This might explain why [29] is closer to ζ eff /s, while [28] is closer to ζ/s.
However one should also note that neither of these theoretical calculations include dynamical resonances like SMASH, and any agreement is probably accidental, as the scattering processes are different in the three calculations.
• Interestingly, the PHSD calculation from Ozvenchuk et al. [31] is not far from our ζ/s, which can be explained in part because PHSD also propagates resonances, and thus including mass changing processes. Using their discrete test particle representation, the bulk viscosity is computed from a discretized version of the relaxation time approximation, where the sum is taken over all particles in the system, which also includes all the resonances.
While this calculation does not account for the dynamical effects of resonances, their widths are explicitly incorporated in the bulk viscosity. In that sense, the effect of long-lived resonances which potentially block the bulk relaxation are also included in the PHSD bulk viscosity calculation.
• The Moroz calculation [32] uses the relaxation time approximation to analytically calculate the viscosities of the hadron gas in a similar fashion as to what was presented with the Chapman-Enskog formalism in Sec. III. In this framework, although all resonances are incorporated in the various cross-sections of the collision term, they do not per se exist as propagating particles in the calculation, and only binary elastic collisions are considered.
As this calculation is closer to our ζ eff /s, we conjecture that the slow processes dominating ζ/s in SMASH are not included in the list of processes of [32], or that the difference is due to dynamical effects not being included in that calculation.
• Finally let us comment on the state-of-the-art values of ζ/s(T ) extracted from hybrid models [68]. The temperature dependence follows some predefined ansatz, motivated by the Hagedorn picture of [23] where ζ/s increases with temperature. A Bayesian analysis is then employed to constrain the functional dependence using experimental data for bulk observables at RHIC and LHC energies. The values of ζ/s extracted at temperatures close to T c come from the 60% confidence intervals in [68] and show a slight increase with temperature. These values are of the same order as our ζ eff /s (which, we remind the reader, should at these high temperatures be considered a minimum value, since some contribution from the higher modes might be missing) but not compatible with ζ/s. This is nicely consistent with the claim that in heavy-ion collisions, the slowest processes (whose inverse rates are larger than the fireball lifetime) do not play any role in the inferred bulk viscosity.

VI. CONCLUSIONS
We have presented our estimates for the bulk viscosity and ζ/s of a hadron gas as a function of temperature between the range T = 80 − 170 MeV at vanishing baryochemical potential. The results at the highest temperatures should be understood as a theoretical extrapolation, as the effects of a deconfined medium-which should take place at such temperatures-are not included.
The calculation of the bulk viscosity is numerically very challenging due to the small size of the fluctuations in the bulk channel, and because the statistical uncertainties in the pressure average can be of the same order. The systematic uncertainty is estimated by comparing to Chapman-Enskog calculations in a simple system with only one particle species. For the final results in a full hadron gas we can confirm that our calculation lies within the same area as previous calculations and extractions from experimental data in heavy-ion collisions. We observe a decreasing trend of ζ/s as a function of temperature, which needs to be reconciled with the expectation of a smooth maximum around the crossover transition to the quark-gluon plasma, which is absent in our model.
We find that mass-changing processes, namely resonance excitations, have a very strong influence on the bulk viscosity. This is rather straightforward to understand since such processes allow to store kinetic energy in the mass of the particles and enhance the fluctuations of the kinetic energy of the system. Our results can be employed in future assumptions for the prior for Bayesian multi-parameter analysis and compared to lattice-QCD calculations once they become available.
One of our main results is the need for a distinction between the inclusive bulk viscosity ζ and the effective bulk viscosity ζ eff . The first one is computed for long-lived systems in equilibrium, in which all components of the medium need to relax for the restoration to equilibrium to occur. We have explicitly shown that the slowest processes determine the bulk viscosity, as their contribution dominates the decay of the correlation function. These modes with relaxation times of several dozens and even hundreds of fm/c make the ζ/s a large coefficient for all temperatures.
The effective bulk viscosity is the coefficient controlling the relaxation to equilibrium of systems with a finite lifetime, as the ones happening in RHICs. The very slow modes then do not have enough time to be relevant and their contribution to the correlation function are explicitly removed. Unsurprisingly, the extraction of ζ/s performed in standard hydrodynamic codes are closer to this effective bulk viscosity, as they match the experimental observables of a finitely lived system.
The effect of these very slow modes was not observed in other coefficients like the shear viscosity, electrical conductivity or cross-conductivities [50][51][52], where the single exponential decay was found to be a good approximation (except for very high temperatures where the system becomes dense). This situation illustrates that the bulk viscosity is a much more subtle quantity than other transport coefficients; being extremely dependent on the microscopical details of the interactions, any comparison between calculations must be performed with caution. The particle and entropy densities can be expressed as n = J 2,1 /T , s = (J 3,1 − µJ 2,1 )/T 2 . (A3) For the quantities used in Eq. (A1) the relations become much more complicated. To simplify the expressions, let us consider the case µ = 0 from now on, which is the one taken in this work. We obtain The adiabatic speed of sound (4) reads In a similar fashion C ζ (0) can be expressed in terms of J n,k (T, µ) functions if desired. We plot some of these thermodynamic quantities as functions of T for different systems containing several hadron species. We consider π, K, N, ρ, K * , ∆, where for the resonances we need to generalize the expression (A2) to include an additional integral over their spectral functions.
In Fig. 12 we present the quantities (A4), (A5) and the adiabatic speed of sound (A6), for a hadron gas when several species are subsequently introduced. Notice that v 2 S already presents a nonmonotonous behavior when ρ mesons are introduced in the pion gas. For a more realistic case with more states covering higher masses, we refer to Fig. 10.
Let us finally comment on two particular cases which can be quite illustrative, although they are not used in the results of this paper. For massless particles, we note that and one obtains ∂P ∂n = 0, ∂P ∂ n = ∂P ∂ S = 1/3, so the bulk viscosity is seen to vanish proportionally to the square of 1/3 − v 2 S [5,66]. For an ensemble where the particle number n is not conserved, one does not introduce any chemical potential, and the thermodynamic functions only depends on T . The speed of sound reduces to, The C ζ (0) in this particular case would read which can be further simplified to Combining the expression for the bulk viscosity in Eq. (7) and our previous result on the shear viscosity [50] (also using the exponential decay ansatz) we obtain where is the shear correlation function at t = 0, and τ η is the relaxation time of a fluctuation in the shear channel. If we assume that τ η τ ζ and introduce the result (A10) for massless particles (E p = p) one gets S 2 p 2 f eq (p) 1 15V d 3 p (2π) 3 p 2 f eq (p) = 15 which coincides with the well-known relation (20) also in the numerical factor. We fit the correlation functions given in Fig. 5 to the form using different methods. First of all, to check that all modes are indeed exponential, we proceed with a sequential method, as described at the end of Sec. IV: one finds the exponential fit of the tail of C ζ (t) and then subtracts the fitted component from the full correlation function. Then, one repeats the procedure to find the exponential decay of the intermediate range of times, and after another subtraction, one fits the small-t part of the function. We present an example of such a fit in Fig. 13 for the temperature of T = 86 MeV (the one with largest error bars). It is difficult to assign an uncertainty to the sequential fit itself, due to the rather manual procedure, so it is given as is.
The quality of the fit is very good. We double-check the resulting fit procedure against a global fit of C ζ (t) using the NonlinearModelFit option in Mathematica [69], and also see that using a larger number N of exponentials results in a poorer quality of the fit, as some components have negative amplitudes, which is physically unreasonable.
The parameters of the "sequential fits" for all temperatures are summarized in Table VI in the upper block of data. All fits have been checked against independent fits in Mathematica (not shown here). Sequential fit T (MeV) C ζ,1 (0) (GeV 2 fm −3 ) τ ζ,1 (fm) C ζ,2 (0) (GeV 2 fm −3 ) τ ζ,2 (fm) C ζ,3 (0) (GeV 2 fm −3 ) τ ζ, 3  We apply yet another method by making a global fit using ROOT [70], which takes into account the error band of C ζ (t) and also provides the statistical uncertainties of the fitting parameters. The outcome of these fits is shown in the lower block of data of Table VI. The numbers are more or less consistent with the "sequential fit", although some deviations remain. Notice that the sequential fit carries an additional systematic error (also difficult to extract) coming from the selection of fit ranges, which has to be decided relatively arbitrarily. In every case, we checked that both independent fits in Table VI describe the correlation function really well, and they actually result in a very similar bulk viscosity. In the main text, the global fit by ROOT is used because it provides a measure of its statistical uncertainty.
For completeness we also provide the results for the global fit in the case of the "smaller" boxes for the same full hadron system. They are shown in Table VII only for the case of the fits using ROOT package.