Non-equilibrium attractor in high-temperature QCD plasmas

We establish the existence of a far-from-equilibrium attractor in weakly-coupled gauge theory undergoing one-dimensional Bjorken expansion. We demonstrate that the resulting far-from-equilibrium evolution is insensitive to certain features of the initial condition, including both the initial momentum-space anisotropy and initial occupancy. We find that this insensitivity extends beyond the energy-momentum tensor to the detailed form of the one-particle distribution function. Based on our results, we assess different procedures for reconstructing the full one-particle distribution function from the energy-momentum tensor along the attractor and discuss implications for the freeze-out procedure used in the phenomenological analysis of ultra-relativistic nuclear collisions.

We establish the existence of a far-from-equilibrium attractor in weakly-coupled gauge theory undergoing one-dimensional Bjorken expansion. We demonstrate that the resulting far-fromequilibrium evolution is insensitive to certain features of the initial condition, including both the initial momentum-space anisotropy and initial occupancy. We find that this insensitivity extends beyond the energy-momentum tensor to the detailed form of the one-particle distribution function. Based on our results, we assess different procedures for reconstructing the full one-particle distribution function from the energy-momentum tensor along the attractor and discuss implications for the freeze-out procedure used in the phenomenological analysis of ultra-relativistic nuclear collisions.
Fluid-dynamic description is a powerful tool in the phenomenological analysis of ultra-relativistic nuclear collisions [1][2][3]. In a fluid-dynamic description of the evolution of the collision system, only a small subset of the degrees of freedom are dynamically evolved. These are quantities derived from the energy-momentum tensor T µν , namely local temperatures, velocities and, in the case of viscous hydrodynamics, information about the shear and bulk viscous tensors. The experiments do not, however, measure fluid-dynamic variables but rather distributions of particles that have "frozen out" and freestream to the detectors -the angular and momentum distributions of these particles inform us about the material properties of the fluid created [4]. To convert fluiddynamic fields to particle distributions, a freeze-out procedure has to be applied. While the energy-momentum tensor depends only on the first momentum-integral moments of the distribution function, the particle distributions contain information about all the moments. Therefore, the conversion of the fluid-dynamic information to particle distributions necessarily involves injection of new information in the form of theoretical assumptions. The common procedure is to assume that the distribution function has a near-equilibrium form whose deviations from equilibrium arise from formally small corrections -the shape of the corrections is determined by the response of a linearized collision kernel in some assumed kinetic theory to an infinitesimal strain [5,6].
As the freeze-out procedure strongly affects the phenomenological analysis and conclusions about the matter created in ultra-relativistic heavy-ion collisions, it is of great interest to scrutinize quantitatively how well justified are the theoretical assumptions about the shape of the non-equilibrium distribution functions. The need for such scrutiny becomes increasingly important in the case of small systems, e.g., peripheral nucleus-nucleus collisions, proton-nucleus, and high-multiplicity protonproton collisions where fluid-dynamical description is being applied to situations which most likely remain far from equilibrium throughout their dynamical evolution.
There has been a large body of work quantifying to what extent various different formulations of viscous fluid dynamics are able to reproduce the time evolution of the components of the energy-momentum tensor undergoing expansion in various geometries. In particular, it has been observed that the hydrodynamic constitutive equations that relate the stress tensor to gradients of the flow fields are well satisfied in systems which are still far from equilibrium at least in systems characterized by flow with a large degree of symmetry [7,8] (for systems with less symmetry cf. [9][10][11][12]) -a feat dubbed hydrodynamization without thermalization . It has also been observed that many microscopic models, as well as various formulations of fluid dynamics, exhibit rapid information loss of some details (in particular the initial longitudinal pressure) of the initial condition leading to non-equilibrium attractor behavior. The qualitative similarity of the attractors between different theories has been advocated to extend the applicability of the fluid-dynamic models to far-from-equilibrium regimes where the ordinary justification of fluid dynamics as a near-equilibrium expansion is questionable.
Much less attention has been paid to the validity of the freeze-out procedure far from equilibrium. One of the challenges is that models with trivial momentum dependence, such as kinetic theory in relaxation time approximation, can give only limited information about the validity of the freeze-out procedure, whereas strongly coupled models without quasiparticle structure do not even possess underlying particle distributions. Here, we discuss the reconstruction of the particle distributions from the energy-momentum tensor in the Effective Kinetic Theory (EKT) for weak-coupling quantum chromodynamics (QCD) that becomes leading-order accurate in the limit of high center-of-mass energy collisions [36]. The rich momentum-dependent structure of the EKT collision kernel allows for a non-trivial test of the freezeout procedure in this theoretically clean limit. We follow 0+1d far-from-equilibrium Bjorken flow within this model and compare different moments of the distribution arXiv:2004.05195v1 [hep-ph] 10 Apr 2020 function to those predicted by hydrodynamic freeze-out prescriptions. We find that EKT seems to exhibit qualitatively similar far-from-equilibrium attractor behavior to RTA kinetic theory and Israel-Stewart-type hydrodynamics [13,17,32] both at early (early-time or pullback attractor) and late times (late-time or hydrodynamic attractor). We find that this attractive behavior is not restricted to the components of energy-momentum tensor but extends to other integral moments as well. We observe that the commonly used freeze-out prescriptions reproduce low-order moments of the distribution well at late times, however, they can fail at early times or when considering moments sensitive to momenta much larger than the temperature. We discuss the implications for phenomenological fluid-dynamic modeling of small collision systems such as pp and pA.
Methodology: We make use of a numerical implementation of the effective kinetic theory (EKT) of Refs. [8,36,37]. In parametrically isotropic systems, EKT gives a leading order accurate description (in α s ) of the QCD time evolution of the one-particle distribution function and allows for a numerical realization of the so-called bottom-up thermalization scenario [38]. In practice, we solve the EKT Boltzmann equation for a gluonic plasma undergoing one-dimensional Bjorken expansion with transverse translational symmetry such that the effective Boltzmann equation reads [39] − where f (p) is the gluonic one-particle distribution function (per degree of freedom). The elastic scattering term C 2↔2 and the effective inelastic term C 1↔2 include physics of dynamical screening and Landau-Pomeranchuck-Migdal suppression and, in order to find the form of the collision kernels, self-energy and ladder resummations are required. For details, we refer the reader to [8,36,37].
For the numerical solution of Eq. (1), we discretize n(p) = p 2 f (p) on an optimized momentum-space grid and use Monte Carlo sampling to compute the integrals appearing in the elastic and inelastic collisional kernels. The algorithm used is based on Refs. [37,40] and exactly conserves energy while also exactly accounts for the particle number violation originating from the inelastic contributions to the collisional kernel. Due to the azimuthal symmetry of Bjorken flow, one can discretize momentumspace on an effectively two-dimensional grid-here we use 250 × 2000 points in the p and x = cos θ directions, respectively. The momenta p are distributed on a logarithmic grid in the ranges [0.02, 45] Λ, where Λ is the typical energy scale of the initial condition. In all Figures presented herein, we used 't Hooft coupling λ = N c g 2 = 10 corresponding to a specific shear viscosity ofη = η/s ≈ 0.62 [14,41].
We follow the time evolution of a complete set of integral moments characterizing the momentum dependence of the distribution function [42] where p = |p|. Note that the energy density is given by ε = νM 20 , longitudinal pressure by P L = νM 01 , and number density by n = νM 10 for ν degrees of freedom (ν = 2d A for d A adjoint colors of gluons). The other moments do not have an interpretation in terms of the usual hydrodynamic moments considered in the literature, although the m = 0 modes are simply related to the effective temperatures introduced in [43,44]. In our results, these moments will be scaled by their corresponding equilibrium values with The temperature T here corresponds to the temperature of an equilibrium system with the same energy density, given by T = (30ε/νπ 2 ) 1/4 . The different moments are sensitive to different momentum regions of the distribution function and for future comparisons, we note that, in equilibrium, the typical momentum contributing to a given moment is , and (c), respectively. These simulations have been initialized with either of the two following initial conditions: (1) spheroidally-deformed thermal initial conditions which we will refer to as "RS" initial conditions [45] and (2) non-thermal color-glass-condenssate (CGC) inspired initial conditions [8]. In the first case, the initial gluonic one-particle distribution function is taken to be of the form where −1 < ξ 0 < ∞ encodes the initial momentumanisotropy and Λ 0 is a temperature-like scale which sets the magnitude of the initial average transverse momentum. In the second case, we take for the form of the initial gluonic one-particle distribution This form has been used in several earlier works (see e.g. [8,10,43,44,46]), and is motivated by the saturation framework, where the initial average transverse momentum scaleΛ 0 is related to the saturation scalẽ [47][48][49]. The overall constant A is set by fixing the initial energy density to match an expectation τ 0 0 = 0.358 νQ 3 s /λ from a classical Yang-Mills simulation of Lappi [49].
In both Figures, the dotted and dashed black lines correspond to EKT evolution using RS-and CGC-type initial conditions, respectively. The integral moments are plotted as a function of a rescaled time variable τ /τ R (τ ), which measures the age of the system in units of the instantaneous interaction time τ R (τ ). As the density of the system changes, so does the interaction time scale, which is given by τ R (τ ) = 4πη/T (τ ). Scaling time in this manner guarantees that, as long as the system is described by hydrodynamics close to thermal equilibrium, M 01 will eventually be described by the first-order gradient expansion, M 01 = 1 − (120ζ(5)/π 5 )τ R /τ , at late times [15,27].
This is independent of the microscopic details or specific values of macroscopic hydrodynamic parameters. This fluid-dynamic behavior is seen in Fig. 1 (a) where the evolutions of all the various different initial conditions eventually converge onto a universal curve -the late-time attractor (see also [50,51]). However, as observed also in simpler toy models, this collapse takes place before the system is well described by the hydrodynamic gradient expansion, the first order of which is shown in Figure 1 (a) as an orange dashed line.
While the late-time attractor behavior for the longitudinal pressure has been observed earlier in simplified kinetic theories, the solutions at hand allow us to study to what extent the attractive behavior determines the full overall shape of the distribution function. Our first main finding is shown in panels (b) and (c) of Fig. 1 which display the time-evolution of two higher moments of the distribution function, M 21 and M 33 . We observe that the higher moments collapse to a universal curve on the same timescale as M 01 , demonstrating that the universality extends beyond simple hydrodynamic moments and it is the entire distribution as a function of p that reaches an attractor form. For corresponding results for a large set of moments see the Supplemental Material associated with this paper. The timescale at which the different solutions to Eq. (1) collapse in Figure 1 is roughly τ /τ R ∼ 0.5. While all theories must eventually collapse on a single curve, the time at which individual solutions collapse to the attractor depends on the details of the model. In [32] two qualitatively different patterns were identified. In RTA kinetic theory and in IS hydrodynamics, the decay to the attractor took place by a powerlaw whose scale was set by the initial time τ 0 such that a unique attractor exists at arbitrarily early time and can be found by studying initialisation with decreasing τ 0 . In contrast, in AdS/CFT, the decay to the attractor happens only at the time scale given by the decay of the quasi-normal modes. Figure 2 shows a set of solutions with fixed initial conditions (RS or CGC) with successively decreasing τ 0 . The figure demonstrates that earlier initializations lead to faster decay to the attractor signifying τ 0 -scaling of the decay and the existence of an early-time (or pullback [29]) attractor in EKT. We note that at late times the attractor for both overoccupied (CGC) and the thermal (RS) initial conditions are the same. This implies that upon reaching the attractor, the late-time evolution is not only insensitive to the initial longitudinal pressure of the initial condition but also to the initial occupancy and momentum profile for these physically motivated initial condition types and moderately large couplings (λ = 10). We note that in each run shown in Fig. 2 we observe a transition from purely free-streaming behavior to a collisionally-broadened longitudinal expan- sion related to the early stages of the bottom-up thermalization [38,50] and to early-time "pre-scaling" behavior seen in prior EKT studies [52]. In Fig. 2 (a) we also compare the EKT attractor to other known attractors for P L /P eq L . The solid purple line corresponds to the exact solution for the attractor in kinetic theory in the RTA approximation [27,32,[53][54][55][56], the orange long-dashed line is the first-order gradient expansion result, the blue dot-dashed line corresponds to a formulation of viscous fluid-dynamics used extensively in phenomenological description of heavy-ion collisions, namely second-order viscous hydrodynamics (vHydro) of Denicol et al [57,58], and the red dot-dot-dashed line corresponds to the anisotropic hydrodynamics (aHydro) attractor [23,[59][60][61][62]. For details concerning how the attractors were determined in each case we refer the reader to the Supplemental Material. We observe that while all of the attractors share some qualitatively similar features, the attractors of the different theories agree quantitatively only at τ /τ R 1 after the attractors follow the hydrodynamic gradient expansion. In particular, we emphasize that in vHydro the longitudinal pressure becomes negative at early times unlike in aHydro or EKT. Panels (b) and (c) of Fig. 2 compare two higher-order moments of the RTA and EKT attractors. While the agreement between the two kinetic theories is rather good at late times for n = 0 moments, the agreement becomes rapidly worse for increasing n [63]. This implies that while the |p|-dependence of the collision kernel may be rather well approximated by the simple RTA kernel at these values of coupling λ, the simplified angular structure of the RTA does not fully capture the shape of the longitudinal structure of the distribution function.
While the fluid-dynamic theories do not specify the higher moments of the distribution functions displayed in panels (b) and (c) of Fig. 2, it is a common practice to infer the full shape of the distribution from the shear components of the energy-momentum tensor only. For a given T µν the linearized viscous correction to the one-particle distribution function, δf can be locally computed given an assumption of the collision kernel. Herein, we consider two possible forms for δf . The (i) quadratic ansatz which results from a wide set of models including RTA with momentum-independent relaxation time, momentum diffusion approximation, scalar field theory, and from EKT in the leading-log approximation [6]. Here Π = Π/ = 1/3 − T zz / is shear viscous correction to the longitudinal pressure scaled by the energy density. At full leading-order, the EKT however has more structure; for large p T , the EKT reduces to power law form of the (ii) LPM ansatz This p 1.5 power-law is numerically close to ∝ p 1.38 that is found to describe the full EKT in the relevant momentum range in [6]. Additionally, we consider a simple (iii) aHydro freeze-out ansatz procedure that does not assume linearization around equilibrium. Instead, in the aHydro freeze-out ansatz, one assumes that the nonequilibrium distribution function can be approximated by a spheroidally-deformed Bose-distribution f (p) = f Bose ( p 2 + ξp 2 z /Λ) [45,59,60]. To test this approach, we fix ξ locally such that the energy-momentum tensor of the ansatz matches with that of the full simulation and then make predictions for higher-order moments.
The different moments obtained by the above prescriptions are compared to the EKT attractor solution in Fig. 3 (see the Supplemental Material for additional moments compared over our entire set of runs). At late times τ > 5 τ R , the low-order moments are described within a few percent by all the prescriptions, while some discrepancy remains even at τ ∼ 20 τ R between the quadratic ansatz (i) and our EKT results. The agreement worsens gradually at earlier times and around τ ∼ τ R where the corrections to longitudinal pressure start to be sizable P L /P eq L ∼ 65%, M 11 exhibits an approximately 20% disagreement between EKT and both linearized ansatze. The disagreement increases for higher moments and for earlier times. In contrast, we observe rather good agreement between the aHydro ansatz and our EKT results at all times.
Conclusions and Discussions: An important step in the phenomenological analysis of nuclear collisions is the freeze-out procedure in which the hydrodynamical fields are converted to particle distributions. In the current phenomenological practice, the quadratic ansatz (i) is widely used. This assumes linear deviations from thermal equilibrium, which is in stark contrast to the far-fromequilibrium conditions in which fluid-dynamical modelling is practiced in current phenomenological applications, in particular in modeling of small systems (see e.g. Refs. [64][65][66][67][68][69]). To address whether these linearized procedures remain quantitatively predictive far from equilibrium, we have confronted them with far-from-equilibrium simulations of QCD effective kinetic theory. Our results in Fig. 3 show that at least in this simplified frameworkadmittedly quite far from the realities of phenomenological modeling -the non-linear aHydro freeze-out ansatz performed better in reconstructing moments of the distribution function compared to linearized ansatze in farfrom-equilibrium systems.
To  [62]. We also implicitly assume that the results as function of τ /τ R do not depend onη as seen for the hydrodynamic moments in [10] and assume a value of η ≈ 0.2, which is consistent with the phenomenological extraction of the quantity.
Using this setup and averaging over our full set of runs, the rescaled times τ /τ R = {0.2, 0.5, 1, 2, 5, 10} map to τ {0.32, 0.86, 1.88, 4.23, 14.1, 38.5} fm/c. This suggests that for τ fo 5 fm/c the lowest-order modes, which are sensitive to p ∼ few T , can be well described by both aHydro and linearized freeze-out prescriptions. This implies that for central ultra-relativistic heavy-ion collisions one can have a faithful reproduction of the lowmomentum part of the freeze-out distribution function. However, when considering higher-moments, which are sensitive to higher momenta p nm , or applying earlytime freeze-out for smaller systems such as peripheral nucleus-nucleus collisions and proton-nucleus collision, the aHydro freeze-out ansatz is favored.
In closing, we note that the current study was performed in a very simple setting with one-dimensional Bjorken flow and considering only massless gluonic degrees of freedom. We leave, for the future, extensions to more realistic geometries [19,20,25] and inclusion of quark degrees of freedom [43,44].