Equilibration and freeze-out of an expanding gas in a transport approach in a Friedmann-Robertson-Walker metric

Motivated by a recent finding of an exact solution of the relativistic Boltzmann equation in a Friedmann-Robertson-Walker spacetime, we implement this metric into the newly developed transport approach Simulating Many Accelerated Strongly-interacting Hadrons (SMASH). We study the numerical solution of the transport equation and compare it to this exact solution for massless particles. We also compare a different initial condition, for which the transport equation can be independently solved numerically. Very nice agreement is observed in both cases. Having passed these checks for the SMASH code, we study a gas of massive particles within the same spacetime, where the particle decoupling is forced by the Hubble expansion. In this simple scenario we present an analysis of the freeze-out times, as function of the masses and cross sections of the particles. The results might be of interest for their potential application to relativistic heavy-ion collisions, for the characterization of the freeze-out process in terms of hadron properties.


Introduction
Kinetic theory [1] has been widely used to study the nonequilibrium evolution of fluids and plasmas, not only for ordinary substances but also in the relativistic domain [2]. For sufficiently dilute systems, the Boltzmann equation (BE) describes how the one-particle distribution function f (t, x, k) relaxes towards equilibrium. Under a general spacetime metric, this equation reads [3] k μ ∂ f (t, x, k) ∂x μ where i λμ are the Christoffel symbols and C [ f ] represents the (nonlinear) collision integral [2]. A non-trivial solution of this equation (aside from the equilibrium distribution) is extremely hard to obtain in this general case. Nevertheless, several approximations can be used to simplify the nonlinear structure of C [ f ]. One of the simplest methods is the relaxation time approximation (RTA) [4], which provides a linearized collision term. In addition, perturbative solutions based on the existence of some small parameter (like the Knudsen number) are also possible e.g. the Chapman-Enskog expansion [5,6]. On the other hand, the BE can also be addressed by pure numerical techniques, known as molecular dynamics simulations or Boltzmann-Uehling-Uhlenbeck (BUU) transport. For relativistic heavy-ion collisions (RHICs), where a mixture of relativistic particles is subjected to mutual interactions and mean-field potentials, the system of coupled BEs can be solved by Monte Carlo methods, as in [7][8][9][10][11][12]. These numerical approaches-more suitable for these complicated systems-also introduce systematic uncertainties, originating from algorithmic approximations and truncations. For this reason, the finding of exact non-trivial solutions of Eq. (1), at least in particular scenarios, is important to test different methods and approximations, either semianalytical or purely numerical.
In RHICs the dynamics and geometry of the created fireball provide certain degrees of symmetry, from which simplified models have been proposed [13]. For some of them, exact solutions of the BE have been found under the RTA. For example, in the Bjorken model [14] (describing a boost-invariant longitudinal expansion) a semianalytic result has already been calculated by Baym [15]. An exact solution in a Gubser expansion (allowing an additional expansion in the transverse plane) has been obtained in [16]. An exact solution in a 3D conformal expanding medium, or Hubble flow, was recently found in [17,18]  This last scenario represents an interesting case with applications to RHICs but also in a cosmological context to describe an expanding universe [19][20][21]. An exact analytical solution of the BE for an expanding medium has recently been presented in [22,23]. This solution is valid for massless particles, interacting in a flat universe under a Friedmann-Robertson-Walker (FRW) metric. The particle-particle interactions are assumed to follow constant total cross sections, but the full nonlinear structure of the collision operator is kept, i.e. no RTA is assumed. The exact solution obtained in [22,23] was used to test a linear approximation of the BE, and approximate solutions based on the RTA.
In this paper we compare a numerical solution of Eq. (1) with this exact result. We employ the new hadronic transport approach SMASH (Simulating Many Accelerated Strongly-interacting Hadrons) [12]. It is used to simulate hot and dense strongly interacting matter with the goal of exploring the quark gluon plasma phase diagram by performing comparisons to experimental data from heavy ion experiments at different accelerators such as the Large Hadron Collider (LHC), the Relativistic Heavy Ion Collider and the SIS-18 at the GSI Helmholtzzentrum für Schwerionenforschung. SMASH constitutes an effective numerical solution of the equations of motion associated with Eq. (1) using a geometrical collision criterion as in UrQMD 1 under a Minkowski metric, i.e. ds 2 = dt 2 − dx 2 − dy 2 − dz 2 . We adapt the SMASH dynamics to work on an expanding homogeneous, isotropic gas of massless particles with a FRW metric ds 2 = dt 2 − a 2 (t)(dx 2 + dy 2 + dz 2 ). 2 In this case the BE is reduced to [21,3,23] where the gain and loss terms present their full nonlinear structure.
Our first goal is to present a non-trivial test of the SMASH code in an expanding geometry by comparing our outcome to the exact solution given in [22,23]. We also check our results against the numerical solution of Eq. (2) for a different initial condition given in [23]. Then, we exploit the flexibility of SMASH to solve the transport equation in an expanding system of massive particles, generating a dynamical freeze-out (or decoupling) due to the Hubble expansion. This opens up the possibility to study more realistic systems of interest in cosmological scenarios, or in RHICs.
In Sec. 2 we present the SMASH solution to the Boltzmann equation for massless particles using several initial conditions, in particular the one for which an exact analytical solution is known. In Sec. 3 we introduce a toy model of freeze-out for relativistic particles when the Hubble rate exceeds the interaction rate. We discuss how the freeze-out time can be extracted from the final spectrum of particles. Finally, we present our conclusions and outlook in Sec. 4.

SMASH solution of the Boltzmann equation under a FRW spacetime
The authors of Ref. [22,23] have calculated an exact solution of the Boltzmann equation (2) for an infinite gas of massless particles with constant elastic cross-section. This is a very particular system, with symmetry properties that help to simplify the transport equation. This scenario is physically motivated by the expansion of the universe in the radiation-dominated era [19].
SMASH is a recently developed transport approach used to describe the hadronic stage of heavy-ion collisions at low and intermediate energies with applications from GSI to LHC physics. In particular, one can use SMASH to simulate a gas of massless particles with constant elastic cross-section σ . In this work we use a spherical volume filled with N particles (this number remains constant due to the absence of number-changing processes). The particles are initialized with an isotropic, homogeneous spatial distribution according to a given initial condition.
Whilst it appears a formidable task to adapt SMASH to a general spacetime, the FRW metric we wish to implement is fairly simple. As SMASH operates in physical phase-space variables, we will always work with the physical 3-momenta k = k phys , as opposed to the approach in [22,23], which uses the covariant momenta. 3 The equations of motion of the particles reflect the physical expansion of the universe. The velocity measured by a comoving observer to the expanding spacetime (sometimes called peculiar velocity) is combined with the Hubble flow: v i is the Hubble parameter. The physical momentum of the particle suffers a redshift and scales as 1/a(t).
Particle collisions are not affected by the Hubble expansion, because the characteristic collision time is always much smaller than H −1 (t), so during the collision the particles do not feel the expansion of the universe.
For a gas of massless particles (radiation) the cosmic scale factor is fixed by the Friedmann equation to be of the form a(t) ∼ t 1 2 [19,21,20]. Following [22,23] we adopt the solution a(t) = (σ denotes the cross section and n 0 the initial particle density) and b r is a parameter which contains the density fraction of radiation in the universe and the Hubble parameter itself at t = 0.
We initialize the particles in a far-from-equilibrium configuration, according to the momentum distribution in [22,23] where λ = exp(μ 0 /T 0 ) is the fugacity of the system, and T 0 a parameter, which can be thought of as an initial temperature of the system. 4 The analytical solution for this initial condition was determined to be where κ(τ ) = 1 − exp − τ 6 /4, and the transformed time variable The distribution function is normalized such that where V is the volume of the sphere, and d 3 k = dkk 2 d k . We show all particles contained in the sub-sphere r < r sub in red, and the particle outside this volume (r > r sub ) in blue (and re-scaled by a factor 1.5 for clarity). We use a single event N ev = 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) We now present the results of SMASH for the solution of the Boltzmann equation for different values of the physical time. Notice that we show these results in terms of k multiplied by the scale factor.
In our simulations we use a constant particle number N = 1.5 × 10 5 contained in a spherical volume with initial radius r 0 = 50 fm. The ratio b r /l 0 is fixed to 0.1. The initial temperature of the system is set to T 0 = 0.2 GeV. With this set of parameters, the fugacity takes a value of λ 0 2.7, which should be constant in the evolution, cf. Eq. (4). To avoid unwanted boundary effects and ensure homogeneity we always work with particles contained in a sub-sphere of radius r sub (t) = a(t)r 0 /1.5. To check the spatial isotropy we plot in Fig. 1 the Mollweide projection of the distribution of particles for t = 0.1 fm. The Mollweide projection reflects the 2D angular distribution of particles (neglecting their distance to the origin) into a 2D ellipsoidal plot. It is widely used in cartography and cosmology, but it can also be applied to heavyion physics [24]. 5 In red we show the particles contained in the sub-sphere r < r sub . In blue we depict the particles with r > r sub , where we have re-scaled the Mollweide projection by the factor 1.5. Isotropy is manifest in the plot.
Before presenting our results for the distribution function we show the time evolution of the particle and energy densities in Fig. 2. The theoretical expressions of the (non-equilibrium) densities admit an analytic form [23], Our results are in very good accordance with these formulas, even for a single event.
In Fig. 3 we present the ratio is the (Boltzmann) equilibrium distribution corresponding to lim τ →∞ f (t, k). 6 We plot the results as a function of ka The output from SMASH is compared to the analytical solution (4) which uses the initial momentum distribution in Eq. (3). The results demonstrate a very good agreement between SMASH and the 5 We thank Felipe Llanes-Estrada for discussions on this topic. 6 While the exact Boltzmann solution cannot be reached at t → ∞ due to a finite τ horizon [23], we still use it for normalization. In Fig. 4 a similar comparison is made with a different initial momentum distribution. We adapt the "two-mode initial condition" introduced in [23]: where L (β) i! are Laguerre polynomials [2]. While no analytical solution for this initial condition has been achieved, one can still numerically solve the BE for it. The latter, when written as a coupled system of equations for the moments of the distribution, can be easily solved by a finite difference method. From Fig. 4 one can see that SMASH results and the numerical solution compare very well. Errors stemming from the binning procedure are more significant here and cause the low-k data points to noticeably deviate from the expected value. Overall this figure, together with Fig. 3, can be considered a validation check for the SMASH approach to solving the relativistic Boltzmann equation in a FRW metric.

Dynamical decoupling of a relativistic massive gas
In this section we present a simple realization for the freezeout of a particle; exploiting the Hubble expansion as a mechanism to dynamically decouple the particles at a certain time. The fundamental idea relies on the fact that the scattering rate |v|nσ 7 decreases in time due to the volume expansion, so that there might exist a time at which this rate becomes smaller than the Hubble rate H and the particles decouple. This is the basic picture behind the description of the thermal history of the primitive universe [19][20][21].
We implement a much simpler scenario, because we only consider a single species of relativistic particles interacting by a constant cross section. The initial temperature is set to T 0 = 0.4 GeV and the mass of the particles to m = 5T 0 , so we deal with a relativistic gas of particles. One should note that a gas of ultrarelativistic or nonrelativistic particles is not useful in this study, as we will see that it is not possible to distinguish an equilibrium distribution evolving in time from a decoupled distribution with a momentum redshift [19,20].
For relativistic particles, it is known that there exists no equilibrium distribution in the collisionless case (e.g. after decoupling). 8 Therefore, our toy model is based on the following scenario: We start with an equilibrium distribution of massive particles. As long as H , collisions allow the system to maintain equilibrium despite of the Hubble expansion. Once ∼ H , the number of collisions will not be enough to re-equilibrate the system, and the distribution function of particles freezes. For later times, when < H the distribution of particles is necessarily out of equilibrium due to the momentum redshift. To be able to see the decoupling within a reasonable amount of time we take the Hubble rate as a pure parameter of the model and use H = βt, with β = 0.01 fm −2 .
Notice that this is a rather arbitrary choice. Our goal, however, is to establish a toy model for the freeze-out in RHICs, with a Hubble parameter encoding the physics of the expanding medium. For this exploratory study, we choose this parameter to illustrate the effect of decoupling in a clear manner.
As said, our initial state is a gas in equilibrium, where the standard Boltzmann distribution of thermodynamics can be applied. The distribution function contains two parameters: the temperature and chemical potential (or fugacity). The latter is needed due to the absence of number-changing processes. Taking into account that collisions maintain equilibrium at early times (when H ), it is possible to compute the values of these two parameters by noting that the evolution of the particle density is 7 |v| is the relative velocity of the colliding particles. 8 More technically, this statement is closely related to the nonexistence of a timelike Killing vector for a FRW metric with a nonconstant a(t). n(T , μ) = n 0 (T 0 , μ 0 )(a 0 /a) 3 , and assuming an adiabatic evolution where the entropy per particle (s/n) is constant in the evolution. These two conditions provide a system of two coupled equations for the values of T and μ. The evolution of the particle density gives T e μ(T ) and the condition s/n = s 0 /n 0 provides μ(T ) where standard thermodynamical relations have been used (see e.g. [25]). It is illustrative to take the ultrarelativistic and nonrelativistic limits of these equations to recover the textbook cases, for m T , (11) which tells us that in an expanding gas of ultrarelativistic (nonrelativistic) particles, the equilibrium state is maintained with a temperature scaling as T ∼ 1/a (T ∼ 1/a 2 ). It is possible to check that we also reproduce the evolution of the chemical potential in these two limiting cases [21].
We check the expected equilibrium evolution using SMASH with relativistic particles with a mass of m = 2 GeV. The cross section is assumed to be constant σ = 40 mb, and the initial configuration is a Boltzmann distribution with temperature T 0 = 0.4 GeV and vanishing chemical potential μ 0 = 0. The geometrical conditions are the same as in the last section. We use 20 events to reduce statistical uncertainty.
In Fig. 5 we present SMASH results at different times together with the expected equilibrium distribution. In the upper panels we check that the SMASH results compare perfectly well to the predicted equilibrium distribution with parameters given by Eqs. (9), (10). In the lower panels we observe that deviations occur, which increase with time. This is the effect of the decoupling, which roughly happens between t = 5 fm and t = 10 fm by looking at the panels. Using the naive condition ∼ H , we numerically obtain an approximate guess of t D ∼ 8 fm. 9 Although this is a very crude estimate, it is consistent with the situation shown in Fig. 5.
In fact, it is possible to determine the form of the distribution function after decoupling. One simply needs to take into account the unavoidable momentum redshift. For any time after decoupling Before fitting our results to this form let us consider the ultrarelativistic and the nonrelativistic limits: These functions mimic the equilibrium distributions with a temperature evolving in the same way as in Eq. (11). Therefore, we 9 is given in SMASH by counting the number of collisions per unit time. recover the result that for these limiting cases, the evolution after decoupling is indistinguishable from an equilibrium situation. 10 Let us take the distribution of our massive particles at some large time t = 20 fm (this would correspond to the detection or measuring time). After checking that the distribution cannot be described by an equilibrium distribution (a χ 2 test over the fit gives values around 6-7 per degree of freedom), we proceed to fit the nonequilibrium distribution to (12). We obtain a very good fit which is shown in Fig. 6. We should stress that the fitting parameters a D and T D are not really independent, but they are related by Eq. (9). However, as the combined implementation of such a nonlinear constraint is involved, we have opted to assume independent parameters, and test the consistency of the results a posteriori by obtaining the freezeout time from each of them.
From the fit we obtain a decoupling temperature of T D = 0.32 ± 0.02 GeV. Using Eq. (9) and the explicit form of a(t) = exp(βt 2 /2) we obtain a freeze-out time of t D = 5.5 +0.8 −0.9 fm. The fit also gives the value of (a/a D ) 2 = 41 ± 4, corresponding to a freeze-out time 10 In a cosmological scenario this claim is not actually true, as the effective degeneracy factor g(T ) of the distribution function is itself a function of temperature, thus producing a difference between the two cases. of t D = 5.2 +0.8 −0.9 fm, which is quite close to the previous number, so we can provide a final average value of t D = 5.3 ± 0.6 fm for the freeze-out time. Notice that this time is consistent with the situation seen in Fig. 5, although smaller than the crude estimate obtained by simply comparing and H . One has to be careful when interpreting these numbers, because all of them assume a sharp freeze-out process, which is certainly not happening in our case, where a smooth decoupling occurs in time. 11 We perform a parametric study of the freeze-out time by varying the mass of the particle and the interaction cross section. We present our results in Fig. 7. From this figure we observe a slight dependence of the decoupling time on the particle mass (for smaller masses this dependence is more pronounced). This can be explained due to the fact that when the mass is increased the relative speed of the particles in the collision, and the particle density, are reduced. The combined effect is a decrease of , which favours the appearance of the decoupling at earlier times. We also see a rather systematic increase of the decoupling time with the total cross section. This is also expected because for larger cross sections increases, producing a delay in the decoupling of particles. If we had taken the condition = H to define the freeze-out times, we would have obtained systematically larger values than those shown in Fig. 7. A similar condition has been proposed in the context of RHICs trading the Hubble rate by the expansion velocity of the fireball surface [26,27]. A slightly different criterion is used in Ref. [28,29], ξ = H , where ξ is a parameter to be determined. Using the values of t D from Fig. 7, we can easily extract the corresponding ξ . We observe that for lower masses ξ is closer to one, but it quickly decreases to ξ 0.3-0.4 for m > 0.5 GeV. As a function of the cross section, ξ is rather stable around ξ 0.4, increasing up to ξ 0.5 for σ = 20 mbarn.
From this simple toy model for decoupling, we can already observe a rather clear dependence of the freeze-out time as a function of the particle properties. Using a more refined model in- 11 Our exponential parametrization of the scale factor helps to have a rather sharp freeze-out process. However, one cannot avoid still having a certain number of collisions even at times as large as t = 20 fm. cluding many interacting species one might address the freeze-out features of RHICs, and verify or falsify some of the hypotheses used to describe this mechanism. In this letter we content ourselves to present the real possibility of such a model and illustrate how the FRW expansion can provide a basis to address the freeze-out of particles in the SMASH transport approach.

Conclusions and outlook
In this letter we have reported our results on the solutions of the Boltzmann equation for a relativistic gas in a FRW spacetime evolution using the SMASH transport approach. For this metric there exists a particular exact solution of the Boltzmann equation which admits an analytical form [22,23]. We have compared the outcome of SMASH using the same initial condition and found a very good agreement for all times. Given that the initial condition is far from equilibrium, this agreement provides a non-trivial check for the SMASH code. We have also presented the solution of the Boltzmann equation for the "two-mode initial condition" introduced in Ref. [22]. In this case no analytical solution is known, but a numerical solution of the Boltzmann equation-written in terms of a system of coupled equations for the moments of the distribution function-is not difficult to obtain. We also find a very good comparison between both approaches for all times.
Thanks to the versatility of SMASH we have solved the transport equation for more general cases. In particular, we have presented here the solution for massive particles, with a Hubble expansion which produces a decoupling of particles when the interaction rate falls below the expansion rate. As long as H , the evolution maintains the equilibrium distribution of particles, thanks to the collisions among them. In this limit we have seen that the predicted equilibrium distributions agree very well with the SMASH outcome. At decoupling (when ∼ H ) the distribution of particles departs from equilibrium due to the absence of collisions. Taking into account the momentum redshift of the particles, we have extracted the value of the freeze-out time which is consistent with the time at which the particle distribution departs from the equilibrium one. Finally we have presented an analysis of the freeze-out time as a function of the mass of the particle and the interaction cross section.
An interesting extension of this "toy model" for the freeze-out is to enlarge the particle content of the system implementing more realistic interactions according to the physics happening in relativistic heavy-ion collisions. This would allow access, in a more systematic manner, to the freeze-out of different species in a quantitative way and explore the hypothesis of a collective decoupling at fixed temperature versus a sequential freeze-out. With the results presented in this letter in a simplified scenario, we have provided indications that a dependence of the freeze-out time on the particles' properties does exist.