Hydrodynamization in systems with detailed transverse profiles

The observation of fluid-like behavior in nucleus-nucleus, proton-nucleus and high-multiplicity proton-proton collisions motivates systematic studies of how different measurements approach their fluid-dynamic limit. We have developed numerical methods to solve the ultra-relativistic Boltzmann equation for systems of arbitrary size and transverse geometry. Here, we apply these techniques for the first time to the study of azimuthal flow coefficients $v_n$ including non-linear mode-mode coupling and to an initial condition with realistic event-by-event fluctuations. We show how both linear and non-linear response coefficients extracted from $v_n$ develop as a function of opacity from free streaming to perfect fluidity. We note in particular that away from the fluid-dynamic limit, the signal strength of linear and non-linear response coefficients does not reduce uniformly, but that their hierarchy and relative size shows characteristic differences.

Most studies of hydrodynamization profit from simplified set-ups that do not reflect all phenomenological complications but that exhibit general features in great clarity. In particular, most studies of hydrodynamization to date assume exact Bjorken boost invariance, employ conformally symmetric collective dynamics and focus on dimensionally reduced 1 + 1D systems [4,5,6,7,9,10,11,12,14,15,16,17,19] (for studies extending this framework, see [13,18,20,21,22,23,24,25,26,32]). Within this setting, one has reached in recent years a thorough understanding of the off-equilibrium evolution of simple observables in various models. For instance, the asymmetry p T /p L between longitudinal and transverse pressure and the higher longitudinal momentum moments of the stress-energy tensor are known to approach rapidly their universal attractor solution in kinetic theory [33,26,16]. The mathematical structures behind this behaviour continue to be studied in the context of resurgence [9,10,27,28].
The lessons learnt from these 1 + 1D systems are expected to carry over to the phenomenological reality in 3 + 1D. For instance, the early-time dynamics of p T /p L in boost-invariant 3 + 1D systems is known to be governed locally in the transverse plane by an effective 1 + 1D evolution, and the 1 + 1D universal attractor for p T /p L is therefore of relevance for the 3 + 1D dynamics. However, very few observables of phenomenological relevance can be studied in 1 + 1D systems, and some important questions have therefore received little attention so far in the debate of hydrodynamization. One of them is whether all bulk observables hydrodynamize under conditions comparable to those under which p T /p L hydrodynamizes, or whether some classes of observables require systems of longer lifetime, larger spatial extent and/or higher density to approach the values they attain under conditions of almost perfect fluidity. Of particular interest in this context are the conditions for hydrodynamization of the azimuthal momentum anisotropies v n of soft multi-particle production, as these are amongst the most abundant and most precisely measured signatures of collective behavior in AA, pA and pp collisions. Here, we analyze their hydrodynamization in a boost invariant conformally symmetric 3 + 1D kinetic transport theory, whose 1 + 1D variants have been used repeatedly in studies of hydrodynamization.
Up until this point, only the linear response coefficients have been studied in full kinetic theory because of the technical challenges related to solving Boltzmann equations for a distribution functions in complex geometries [37,25,26], though some results exist for perturbative solutions around free-streaming [38,37]. We have developed numerical techniques to solve such systems and we present here the first non-linear response coefficients, and we present the first solution to the Boltzmann equation for an initial condition with realistic event-by-event fluctuations.
Kinetic Theory. We consider massless, boost-invariant kinetic theory in the isotropization-time approximation, and we restrict the discussion to the first momentum moments F( x ⊥ , Ω, τ) = 4πp 2 d p (2π) 3 p f of the distribution function f . Here, p is the modulus of the three-momentum, the velocity is v µ ≡ p µ /p with p µ p µ = 0 and v 0 = 1, and Ω denotes the angular phase space of v µ . F defines the energy momentum tensor T µν = dΩ v µ v ν F, as well as arbitrary higher v µ -moments that lie beyond hydrodynamics. It satisfies the equations of motion [26] where ε is the local energy density. Fluid-like and particle-like excitations are known to coexist in this kinetic transport and their properties can be calculated analytically. In particular, the coupling γ is related to the specific shear viscosity η sT = 1 5γε 1/4 , and F relaxes locally on a time scale τ R = 1 γε 1/4 to the isotropic distribution F iso (τ, x ⊥ ; Ω) = ε(τ, x ⊥ ) (−u µ v µ ) 4 whose functional form is fixed by symmetries and by the Landau matching condition, u µ T ν µ = −εu ν . As the dynamics (1) is scaleless, dimensionful characteristics of the collision system can enter only via the initial conditions, and they can affect results only in dimensionless combinations. For a system of transverse r.m.s. size R and energy density ε 0 at initial time τ 0 , it follows that the opacityγ = γR 3/4 (ε 0 τ 0 ) 1/4 is the unique model parameter. Eq. (1) interpolates between free-streaming in the limit of vanishing opacityγ → 0 and ideal fluid dynamics in the limitγ → ∞.
We initialize (1) with two different classes of initial conditions. We first study linear and non-linear response coefficients based on the simple Gaussian ansatz The exponential exp − r 2 2R 2 multiplying the cos-term ensures that the distribution stays positive everywhere for sufficiently small δ n 's. The initial spatial azimuthal asymmetries are proportional to the real factors δ n , and they are oriented along the azimuthal directions ψ n . Alternatively, we initialize (1) also with the "realistic" initial conditions arising from the T R ENTo model by replacing the radial profile with that arising from the initial state model.
For both classes of initial conditions, we quantify azimuthal anisotropies in terms of the complex-valued spatial eccentricities for n > 1, Evolving with eq. (1) the initial conditions (2), we obtain the evolution of the energy-momentum tensor T µν and the transverse energy flow dE ⊥ at late times This determines the energy flow coefficients V n = v n e i n φ n , where φ n is the azimuthal orientation of the energy flow. In contrast to flow coefficients extracted from particle distributions dN, our study focusses on energy-flow coefficients which are not affected by hadronization since hadronization conserves energy and momentum. The viscous fluid-dynamic limit of eq. (1) is restricted to the evolution of seven fluid-dynamic fields which may be identified with those seven components of T µν (τ, x ⊥ ) = dΩ v µ v ν F that do not vanish under boost-invariance. We are interested in the apparently simple kinetic theory (1) for F(τ, x ⊥ ; φ, v z ) away from the fluid dynamics limit since it provides an explicit realization of fluid fields coupled to a tower of arbitrarily many non-fluid-dynamic excitations (that may be parametrized by the higher v µ -moments of F). However, going beyond the fluid-dynamic limit has a price: F depends on two additional dimensions φ and v z in momentum space. Discretizing φ in twenty points and discretizing the v z -dependence in 50 points implies a 1000-fold increase of the numerical complexity compared to viscous fluid dynamics. The numerical method for solving this evolution equation (1) has been described in [26], but there it was applied only to the linear response of flow coefficients for infinitesimally small n when the coupling between different harmonics can be neglected and the numerics simplifies. Here, we overcome this remaining limitations and we study the kinetic theory for arbitrary eccentricities, arbitrary opacities, and arbitrary coupled non-linear responses.
Results for mode-by-mode kinetic theory. The coefficients v n are known to arise from the dynamical response to spatial eccentricities n in the inital nuclear overlap. The numerically largest responses are linear (v n ∝ w n;n n ) [34], but sizeable quadratic (∝ w n;n 1 ,n 2 n 1 n 2 ) and cubic corrections have been quantified [35,36] and these can dominate higher harmonics (n ≥ 4). For linear responses to spatial eccentricities, there is an intrinsic ambiguity between the initial geometry that specifies the values n , and the collective dynamics that builds up v n from these n . Non-linear response coefficients are of particular interest, since they help to disentangle this ambiguity.
As a first example, we consider initial conditions (2) in which a single mode δ 2 is excited (δ n = 0 for n 2). In the course of the evolution, the non-linear mode-mode coupling of this initial second harmonic with itself excites 4th, 6th, 8th, ... harmonics, but also the 0th harmonic. In turn, these higher harmonics affect the non-linear response of v 2 . For this reason, numerical studies of the non-linear response to v 2 require sufficiently fine discretization in the momentum angle φ to follow numerically also the higher excited harmonics. The numerical results shown here were obtained for a φ-range discretized with 40 points, and their numerical stability was checked with finer discretizations. Our first main result is to observe that the non-linearities are more important for large opacity, as the lines in Fig. 1 develop larger slopes and curvatures. While the numerical results for v 2 and v 4 in Fig. 1 do not involve a perturbative expansion in n orγ, symmetry arguments imply that they must agree for sufficiently small 2 with the perturbative series V 2 = w 2,2 2 + w 2,222 2 * According to eq. (5), the response coefficient w 2,2 (γ) at a given opacityγ is the intercept lim 2 2 →0 v 2 2 ( 2 2 ) of the corresponding curve in Fig. 1 with the ordinate. The non-linear response coefficicient w 2,222 (γ) is the slope of the same curves in Fig. 1  ( 2 2 ). For notational simplicity, we do not denote explicitly the phases of the eccentricities in the following as these can be inferred easily from symmetry arguments. Fig. 2 shows theγ-dependence of the linear response coefficient extracted from Fig. 1 in this way.
In close analogy, we determine other linear and non-linear response coefficients numerically by seeding the initial We ask next how the linear and non-linear response coefficients in Figs. 2 and 3 hydrodynamize, i.e., how they approach their fluid-dynamic limit with increasing opacityγ. To this end, we relate the opacity that characterizes kinetic transport to quantities accessible in viscous fluid dynamics. The definitionγ = γ R 3/4 (ε 0 τ 0 ) 1/4 assumes that the early-time evolution is given by free-streaming which is not the case for viscous fluid dynamics. We therefore have to work with an equivalent definition that can be expressed in terms of quantities measured at a time at which the flow builds up and fluid dynamics may be operational. To this end, we writê where, for the Gaussian background in the initial condition (2), ε 0 and ε R denote central (r = 0) energy densities at times τ 0 and R, respectively. The function f 0→R (γ) = ε R R ε 0 τ 0 is defined as the ratio of the energy per unit rapidity at time τ = R to the energy which the system would have if it were free-streaming [26]. We calculate f 0→R (γ) from kinetic theory forγ ≤ 10, and we match for largerγ to the known asymptotic large-γ behavior f 0→R ∼γ −4/9 .
With f 0→R (γ) known, we relate viscous fluid-dynamic calculations toγ by specifying ε R and η/s from fluid dynamics and solving eq. (6) forγ. In particular, we use the kinetic relation between the interaction strength γ and the shear viscosity η sT = 1 5γε 1/4 . We then initialize at some initial time τ 0 < R the components of T µν from the same initial conditions (2) as the kinetic theory and we evolve them with viscous fluid dynamics for varying η s and ε 0 . This allows us to determine ε R andγ, and to extract from the transverse energy flow (4) at late times the energy-flow coefficients v n . In general, these results depend on τ 0 . That the τ 0 → 0-limit of v n (γ) exists is a direct consequence of the fact that viscous fluid dynamics, like kinetic theory, has a universal attractor solution at arbitrarily early times [12]. While the attractor of kinetic theory keeps ετ fixed leading to the scaling ofγ = γ R 3/4 (ε 0 τ 0 ) 1/4 , the attractor of the viscous (Israel-Stewart) hydrodynamics considered here keeps ετ 4 15 √ 5−5 constant [13]. Therefore taking the τ 0 → 0 limit while keepingγ as defined in eq. (6) fixed corresponds to scaling initial energy densities by ε(τ 0 ) ∝ τ (which differs between kinetic theory and the fluid dynamics) for the physical observables measured in experiments and it informs us about the extent to which the entire signal v n is or is not build up by the degrees of freedom encoded in viscous fluid dynamics. For linear response coefficients, this comparison is shown in the right panel of Fig. 2. Technically, we evolve the viscous fluid-dynamic equations as described in Ref. [39,40] by splitting all fluid dynamic fields into an azimuthally symmetric background and an azimuthally anisotropic perturbation and solving for them to first order in initial eccentricites. In the same way, we set up a control calculation for the much simpler ideal fluid-dynamic equations to obtain an independent determination of linear response coefficients in the limitγ → ∞ (arrows in Fig. 2). Results for w 2,2 and w 3,3 differ somewhat from those reported in [26] since the initial conditions are different.
As expected from general reasoning, the viscous fluid-dynamic results for v n (γ) n in the limit τ 0 → 0 asymptote forγ → ∞ to the ideal fluid-dynamic results in the same τ 0 → 0 limit, see Fig. 2. Remarkably, the hierarchy between the elliptic and triangular linear response coefficient gets inverted as a function ofγ: kinetic theory at lowγ shows w 2,2 > w 3,3 while ideal fluid dynamics shows w 2,2 < w 3,3 . Viscous fluid dynamics accounts for this inversion qualitatively: for very small specific shear viscosity η s , i.e., very large opacityγ, it is consistent with ideal fluid dynamics, but the hierarchy changes as a function of opacity, see right panel of Fig. 2. Also the results from kinetic theory hint at such an inversion, as the slope of w 3,3 (γ = 10) is larger than the slope of w 2,2 (γ = 10).
As seen from Fig. 2, viscous fluid dynamics reproduces the main qualitative trends of kinetic theory (hierarchy of response coefficients) atγ ∼ O(10), but significant quantitative differences persist. On general grounds, we expect that kinetic theory matches quantitatively to viscous fluid dynamics at sufficiently largeγ when the fluid dynamic gradient expansion becomes quantitatively reliable. All data shown here are consistent with this expectation. It would clearly be interesting to extend the numerical calculations in kinetic theory to largerγ and to determine theγ-scale at which a seamless matching to viscous fluid dynamics is found. However, with increasingγ, the numerical evaluation becomes more expensive, and within the scope of the present letter, we were not able to push to higherγ.
We have extended this analysis to a set of quadratic and cubic response coefficients, see Fig. 3. To make some statements about their hydrodynamization we determine the quadratic response coefficients in the limitγ → ∞ by solving ideal fluid dynamics to second order in eccentricities (arrows in Fig. 3). Within the rangeγ < 10, several quadratic response coefficients are seen to cross, and atγ = 10, the hierarchy of the numerically large response coefficients (w 5,23 > w 4,22 > w 6,24 ) found in kinetic theory is consistent with that of ideal fluid dynamics. In the rangê γ > 10, the numerically smaller response coefficients w 3,25 and w 2,53 need to cross. These observations give further support to the conclusions reached from Fig. 2.
In a remarkable note [38], it was observed already that in the dilute limit of kinetic theory far from equilibirum, linear and quadratic response coefficients grow linearly in the average number of rescatteringsN resc. while cubic ones have a quadratic dependence. In Ref. [38], this scaling was established for elastic two-to-two collision kernels. The line of arguments of Ref. [38] does not apply to the collision kernel (1). However, a perturbative expansion of (1) inγ can be viewed as an expansion in the average number of scattering centers [37], and it is therefore natural to test whether our results show this same scaling, too. For linear and quadratic coefficients, we know already from the perturbative analysis in [37] that they do. For cubic response coefficients, however, we observe small violations of the scaling. In the neighborhood ofγ = 0, the cubic coefficients in the right panel of Fig. 3 show a small linear component, though the quadratic one can be dominant.
Evolving initial conditions with realistic event-by-event fluctuations in kinetic theory. We now apply our newly developed numerical machinery to the first exploratory study of a realistic initial condition that would be one single event in an event sample of an event-by-event analysis. The initial condition is a typical T R ENTo event [41] in the 5 − 10% centrality class smoothened such that only initial n 's for n ≤ 7 are kept. We have checked that the finesse of our discretization allows for the stable propagation of such events. A typical time evolution is shown in the upper panel of Fig. 4 withγ = 2. It illustrates that the Boltzmann equation can be solved non-perturbatively for distribution functions representing realistic initial conditions.
The radial profile of the T R ENTo event studied here differs from (2) and this can affect the value of linear and non-linear response coefficients. To quantify the difference, we compare the w 2,2 extracted for these two profiles and find the following numbers w (TrENTo) is not a linear response coefficient, since it was extracted at finite eccentricity, but Fig. 1 informs us that the numerical contribution arising from finite eccentricity is negligible for small opacity. We checked this for the T R ENTo profile as well (data not shown). We observe that the dependence on the radial profile in the linear response coefficients ranges from 5% to 15% in thiŝ γ-range. The analogous study of w 3,3 shows a 2% to 10% difference in the sameγ-range. Therefore, the open circles in the lower panel of Fig. 4 are accounted for within 2% -15% accuracy by the linear response coefficients calculated from the simplified profile (2). The remaining difference between open circles and full results in Fig. 4 result from mode-mode couplings of different harmonics. We see that while the linear response covers the ballpark of the results, non-linearities have to be included to go reliably beyond 20%-30% accuracy. The non-linearities generated by the lowest harmonics n ≤ 3 account for half of all the non-linearities. This paper is motivated by the wealth of studies of hydrodynamization and thermalization in simplified settings. We have developed the necessary machinery for overcoming many of these simplications and to facilitate studies of hydrodynamization in complex realistic geometries, and to thus push the study of hydrodynamization from in vitro to in vivo. The ability to solve the Bolzmann equation for ultra-relativistic systems with realistic initial geometries and including all non-linear mode-mode couplings provides insight into how the characteristic features of fluid dynamics emerge gradually with increasing interaction strength. Away from the fluid dynamic limit, signals of collectivity are not simply reduced uniformily in size, but their relative strength varies characteristically with opacity, the hierarchy of the dominant linear response coefficients is inverted and so is the hierarchy of several non-linear ones. This may provide novel possibilities for characterizing to what extent systems of different size do or do not hydrodynamize. In the long run, we hope that the technical advances documented here can be developed further to study the evolution of event samples, and to study Boltzmann equations with other phenomenologically relevant complications.