Flow in AA and pA as an interplay of fluid-like and non-fluid like excitations

To study the microscopic structure of quark–gluon plasma, data from hadronic collisions must be confronted with models that go beyond fluid dynamics. Here, we study a simple kinetic theory model that encompasses fluid dynamics but contains also particle-like excitations in a boost invariant setting with no symmetries in the transverse plane and with large initial momentum asymmetries. We determine the relative weight of fluid dynamical and particle like excitations as a function of system size and energy density by comparing kinetic transport to results from the 0th, 1st and 2nd order gradient expansion of viscous fluid dynamics. We then confront this kinetic theory with data on azimuthal flow coefficients over a wide centrality range in PbPb collisions at the LHC, in AuAu collisions at RHIC, and in pPb collisions at the LHC. Evidence is presented that non-hydrodynamic excitations make the dominant contribution to collective flow signals in pPb collisions at the LHC and contribute significantly to flow in peripheral nucleus–nucleus collisions, while fluid-like excitations dominate collectivity in central nucleus–nucleus collisions at collider energies.


Introduction
Azimuthal asymmetries v n in soft transverse multi-particle distributions have been measured in detail in PbPb, pPb and pp collisions at the LHC [1][2][3][4][5][6]. At lower energies they have been studied in nucleus-nucleus and deuteron-nucleus collisions at RHIC [7][8][9], and at yet lower center-of-mass energies reached in fixed-target experiments, e.g., [10]. For sufficiently central nucleus-nucleus collisions at collider energies, model analysis support a fluid-dynamic interpretation of the measured asymmetries [11,12]. The recently estaba e-mail: a.k@cern.ch b e-mail: urs.wiedemann@cern.ch c e-mail: b.wu@cern.ch lished persistence of sizeable azimuthal asymmetries in the smaller pPb and pp collision systems [3][4][5][6] now raises the question of whether small transient fluid droplets are also formed in these systems, or whether other physics mechanisms are at the origin of the observed signatures of collectivity. The main purpose of the present work is to address a particular variant of this central question: To what extent are fluid-dynamic or particle-like excitations at the origin of the flow phenomena observed in proton-proton, proton-nucleus and nucleus-nucleus collisions? And how does the interplay between these two sources of collectivity change as a function of system size and energy density?
To set the stage for these questions, we recall that any self-interacting matter that does not spontaneously break Lorentz symmetry carries both fluid-dynamic and non-fluid dynamic excitations. Analysis of the retarded two-point correlation functions G μν,αβ R (x; t) = [T μν (x, t), T αβ (0, 0)] is one way of characterizing the excitations that a system can carry. In particular, the Fourier transform G μν,αβ R (k, ω), understood as a function of complex frequency ω, displays non-analytic structures that can be related unambiguously to excitable physical degrees of freedom. Fluid dynamic excitations correspond to poles at positions ω hyd (k) close to the real axis in the negative imaginary ω half plane. In all known stable and causal dynamics, they are supplemented by other non-analytic structures whose precise nature depends on the microscopic dynamics. Prominent examples are sketched in Fig. 1. For instance, for the non-abelian plasma of N = 4 super-Yang-Mills theory in the large-N c limit at strong coupling, G μν,αβ R (k, ω) displays-in addition to fluid dynamic poles-a tower of quasi-normal modes that are located deep in the negative complex plane at depths ∝ −i2π T n, n ∈ 1, 2, 3, ... [13][14][15]; somewhat distorted versions of this non-fluid dynamic pole structure are found in a class of related, strongly coupled quantum field theories with gravity duals [16][17][18]. In contrast, kinetic theory in the relaxation time approximation (RTA) supplements fluid-dynamic The presence of hydrodynamic poles in the long wavelength limit is a common feature of all Lorentz symmetric theories with selfinteractions. On the contrary, the structure of the non-hydrodynamic sector closely reflects the underlying microscopic degrees of freedom excitations with a branch-cut of quasi-particle excitations that is located at a depth −i/τ R set by the relaxation time τ R [19]; scale-dependent versions of this RTA allow for more general but related non-analytic structures [20]. Also Israel-Stewart dynamics [21] is more than viscous fluid dynamics. The ad-hoc requirement that the shear-viscous tensor relaxes to the constitutive fluid-dynamic 1st-order relation on a timescale τ π results in a non-fluid-dynamic mode in Israel-Stewart dynamics.
Very little is known about the possible excitations in finite-temperature QCD at energy densities attained in heavy-ion collisions. In the limit of vanishing interactions, G μν,αβ R (k, ω) has a tower of branch cuts corresponding to different Matsubara frequencies [15], but one may only speculate whether these structures evolve gradually with the strength of α s .
The experimental determination of these non-hydrodynamic structures would be tantamount to determining the microscopic structure of quark-gluon plasma. While experimental reality prevents us from directly measuring the response of a static and infinite quark-gluon plasma to a linearized perturbation and directly accessing the nonhydrodynamic structures in Fig. 1, the non-hydrodynamic sector may leave its imprint to the dynamical evolution of the system created in physical collisions. This is especially the case if the hydrodynamic and non-hydrodynamic modes are not well separated from each other, or even more so if the non-hydrodynamic modes are closer to the real axis than the hydrodynamic ones and thus dominate the evolution of the system.
The non-hydrodynamic sector becomes more prominent in smaller systems. Firstly, the smallest wavenumber that excitations in a system with transverse size R can have is k ≥ 1/R. With decreasing R, the spectrum of fluiddynamic excitations in the system thus corresponds to poles that lie deeper and deeper in the negative imaginary ω-plane (at positions Im ω hyd − η ε+ p 1 R 2 for the shear channel), and thus approach the non-hydrodynamic sector. Secondly, in a system that lives for a lifetime of τ , the information of modes with Im(ω) > 1/ τ is lost as the decay of the modes is dictated by the imaginary part of frequency exp(Im(ω)τ ). Therefore with decreasing τ ∼ R, one may access an increasing amount of information about the nonhydrodynamic structures that lie deeper in the complex plane. For both of these reasons, systems that are small enough are inevitably dominated by physics that is not fluid-like and the specific way of how the hydrodynamic description fails for small systems carries the information of what lies beyond.
A systematic comparison of system-size dependencies between models of different microscopic dynamics and measured azimuthal anisotropies may offer an empirical pathway to discriminate between different possibilities. We note that while disagreement with the data would exclude a given model, the opposite does not hold true. In fact, an example of such an agreement comes from the phenomenological success of Israel-Stewart hydrodynamics in small systems where the dynamics is dominated by the non-hydrodynamic pole governing the relaxation of the shear-stress tensor to its Navier-Stokes value. Apparent agreement with data surely does not imply that the microscopic dynamics of Israel-Stewart theory are at play in Nature but does mean that the non-hydrodynamic sector of Israel-Stewart theory at least somewhat resembles that of the true underlying dynamics.
To what extent we may discriminate whether the underlying dynamics contains quasiparticles, quasi-normal modes of AdS black hole, or perhaps something else depends to what extent we can identify features inconsistent with data with decreasing system size. The most prominent feature exhibited by the small systems is the free-streaming dynamics upon which the standard modelling of pp-collisions is based. Hence, it is a natural starting point to ask how this limit is approached in the only one of the above mentioned models that does support free streaming, namely the kinetic theory.

Background
Kinetic transport has been considered since the beginning in the study of relativistic nucleus-nucleus collisions [22]. On the one hand, there are parton-cascade formulations [23][24][25][26][27] that follow the time evolution of the individual partons. On the other hand, direct solutions of the Boltzmann transport equation trace the particle distribution functions.
Parton cascades are nowadays embedded in widely used phenomenological models like the AMPT [28] and BAMPS [29] that aim at a complete dynamical description of all stages of nucleus-nucleus collisions, and that reproduce in partic-ular many aspects of the measured azimuthal anisotropies in nucleus-nucleus collisions, and their system size dependence [30][31][32][33]. Parton cascades have also been used to study the Boltzmann transport equation in isolation with the aim of gaining qualitative insights into how kinetic theory approaches a hydrodynamic regime [24,26,[34][35][36], how it builds up collective flow [37][38][39][40], and how this affects the flow of heavy quarks [41]. In addition to the above mentioned partonic cascades, we mention in passing that hadronic transport codes, such as [42][43][44], play an important role in phenomenology.
Aspects of real-time dynamics of QCD can be followed in an effective kinetic theory (EKT) [45] which provides a systematic expansion in the strong coupling constant α s . Transport coefficients of weak-coupling QCD have been evaluated by directly solving the corresponding Boltzmann transport equation near thermal equilibrium [46][47][48][49][50][51]. Solutions to the EKT Boltzmann equation have been studied also in more complex geometries. In particular, solutions to the EKT Boltzmann equation in 1 + 1 dimensional longitudinally expanding systems have shown that the out-of-equilibrium systems of saturation framework (see, e.g., [52])-described by classical gluodynamics at early times [53]-can be smoothly evolved to late times where the system allows for a fluid-dynamic description [54][55][56][57]. In support of its potential phenomenological relevance, this perturbative approach-if supplemented with realistic values of the coupling constantleads to short hydrodynamization timescales [55] comparable to those obtained from non-perturbative strong coupling techniques [58][59][60][61][62] and favored in phenomenological models. This approach has been extended to 2 + 1 dimensional solutions [63] which form the foundation of tools that can be used to match the out-of-equilibrium initial stage models to a hydrodynamic stage in nucleusnucleus collisions [64,65]. In these solutions the transverse translational symmetry is, however, broken only by small linear perturbations limiting their use in the study of small systems which potentially do not reach the fluiddynamic regime before disintegrating due to transverse expansion.
In the past years there has been significant interest also in solutions to Boltzmann transport equations with simplified collision kernels, such as relaxation time approximation or diffusion-type kernels [66][67][68][69]. These formulations have the advantage that due to their simpler structure they allow for more detailed analysis. That has allowed the study of, e.g., the role of out-of-equilibrium attractors and their relation to the non-hydrodynamic modes [66,70]. While this class of models is not directly based on QCD (cf. [51,71]), they do share characteristics with EKT and may be used as toy model of EKT. In addition, they may reveal qualitative features of systems with quasiparticle descrip-tions even in a regime where weak coupling approximations of EKT are not appropriate. Besides works that solve the Boltzmann transport equation perturbatively in the collision kernel around the free-streaming solution [40,72,73], these studies are still constrained to simple geometry of 1 + 1 dimensions with the notable exception of [74] that breaks the translational symmetry but retains a de Sitter symmetry and remains restricted to azimuthal rotational symmetry.
The first non-perturbative solution to Boltzmann transport equation under longitudinal transverse expansion and broken transverse translational symmetry, including broken azimuthal symmetry was given in [75]. The approach in [75] that we document in detail in Sect. 3 one starts from an essentially analytic formulation of kinetic transport in isotropization time approximation (ITA) and extends it to the calculation of azimuthal flow. The fluid-dynamic properties of the model are analytically known which allows us in Sect. 4 to disentangle unambiguously and quantitatively fluid-dynamic from non-fluid dynamic features in the dynamical evolution of kinetic theory. We thus investigate a system in which the modern conceptual debate about a quantitatively controlled understanding of hydrodynamization that has remained largely limited so far to arguably academic, lower dimensional analytic set-ups, can be extended to a problem of sufficient complexity to provide insight into the current phenomenological practice. We hasten to remark that the kinetic transport theory considered here remains a simple one-parameter model designed mainly to address conceptual questions like those posed at the beginning of this introduction. However, in an effort of pushing the relevance of this starting point as far as possible into the realm of full phenomenological studies explored so far with multi-parameter fluid dynamic and parton cascade codes of significant complexity, we shall compare in Sect. 5 results of this one-parameter kinetic theory to results on the system size dependence of measures of collectivity in hadronic collisions at RHIC and at the LHC.

The model
We solve the kinetic theory in isotropization-time approximation (ITA) with boost invariant geometry but with strongly broken transverse translational symmetry. Our initial conditions have an approximate azimuthal symmetry that is broken by linear perturbations. In this section, we describe the model that we solve using free-streaming coordinate-system methods introduced in [75] that we fully document in Appendix Appendix A.

Isotropization-time approximation
Our starting point is the longitudinally boost-invariant massless kinetic transport equation [22] Here the distribution function f (τ, x ⊥ ; p ⊥ , p z ) denotes the phase space density of massless excitations of momentum In the following, we will restrict our discussion to observables that can be obtained from the first p-integrated moment of the distribution function which does not depend anymore on the magnitude of the momentum p, but does depend on the direction of the propagation of particles v ⊥ and v z collectively denoted by , subject to constraint that v 2 z + v 2 ⊥ = 1. The evolution equation of the integrated distribution function is obtained by taking the first integral moment of Eq. (1) This equation closes if the collision kernel C[F] can be expressed in terms of F. While this does not generically happen, it is the case in the isotropization-time approximation (ITA). The ITA is a simple implementation of the physical picture that interactions drive the system towards isotropy even if the system is far from equilibrium and may be too short-lived to thermalize. The ITA therefore assumes that the distribution F evolves in the local rest frame u μ (τ, x ⊥ ) to an isotropic distribution F iso in timescale τ iso , and that the collision kernel is taken to be of the form with v μ = p μ / p = (1, v), and where the local rest frame is given by the Landau matching condition u μ T ν μ = −εu ν . We impose conformal symmetry that guarantees that the isotropization time is related to the local energy density ε The constant of proportionality γ is the only model parameter in the ITA. In particular, the form of the isotropic distribution is completely determined by Lorentz and conformal symmetry.
As discussed already in previous work [73,75], the isotropization time approximation is found to reproduce the evolution of the QCD weak coupling effective kinetic theory within ∼ 15%. It shares important commonalities with the well-known relaxation time approximation, but it is free of the assumption that the distribution function f evolves directly to a thermal equilibrium. As such it is better suited for the description of non-equilibrium systems such as anisotropic plasmas to which a local equilibrium distribution is difficult to associate to.

Initial conditions
We specify initial conditions for the first momentum moments F(τ 0 , x ⊥ ; φ, v z ) of the distribution function at initial proper time τ 0 → 0, with the azimuthal momentum angle φ defined by v ⊥ = (cos φ, sin φ). As particle production is local in space, the initial momentum distribution is azimuthally symmetric at initial time τ 0 for statistical ensembles of collisions. Therefore, F cannot depend on φ initially. Moreover, since ultra-relativistic pp, pA and AA collisions are expected to display very large momentum anisotropies between the transverse and longitudinal direction, the initial condition will be taken to be of the form The initial transverse spatial x ⊥ -distribution is given by an azimuthally isotropic profile function P iso (r/R), with a characteristic root mean square transverse system size R = d 2 r r 2 P iso / d 2 r P iso 1/2 . The isotropic profile is multiplied by an azimuthal dependence with small δ m,n P aniso where θ is the coordinate space azimuthal angle x ⊥ = r (cos θ, sin θ). The azimuthal eccentricity of the initial spatial distribution is characterized by the harmonic coefficients δ m,n that are related to the initial state eccentricities m (cf. Eq. (30)). The initial conditions are then of the form We consider two transverse spatial profile functions P iso : • Initial Gaussian profile The normalization of this isotropic distribution is chosen such that the initial energy density ε(τ 0 , r = 0) calculated from F(τ 0 , x ⊥ ; φ, v z ) in Eq. (9) satisfies ε(τ 0 , r = 0) = ε 0 . • Initial Woods-Saxon (WS) profile For a Pb nucleus, the standard Woods-Saxon parametrization of the density profile is ρ WS (r, z) = 1/ exp where R WS = 1.12 A 1/3 − 0.86A −1/3 | A=208 = 6.49fm, and δ WS = 0.54fm. The corresponding nuclear profile function T WS (r ) ≡ ∞ −∞ ρ WS (r, z) dz is the projection of this density onto the transverse plane. Simple models of the local energy density in the transverse plane take ε(τ 0 , r ) proportional to the nuclear overlap of two Pb nuclei, i.e., for vanishing impact parameter, proportional to T 2 WS (r ) [76]. This motivates us to make the ansatz where we choose the number c rms ≈ 4.42 such that the radius mean square of P iso r R is properly normalized to unity.
In principle, results can depend on the initialization time τ 0 . However, at sufficiently early times, the system must evolve close to free streaming, and varying the initial time will therefore not matter as long as ε 0 τ 0 is kept constant. The τ 0 → 0 singularity of the initial conditions (A56) is in this sense benign since it can be scaled out of the solutions of kinetic theory for constant ε 0 τ 0 (For more detailed discussion see, Appendix A3). This allows us to present results in terms of a single physically meaningful model parameter, the opacityγ As explained in detail in Sect. 5.6, this parameter is related to the transverse size of the system in units of the mean free path, R τ iso . It is linearly related forγ 1 but in opaque systems the relation is R τ iso ∝γ 8/9 .

Numerical results
We analyze now in detail how the time evolution of the energy-momentum tensor varies with transverse system size between the extreme limiting cases of free streaming and ideal fluid dynamics by varying the only available dimensionless parameter, the opacityγ . In the following, at z = 0 time t and radial r coordinates are plotted in units of the system size R.

Time evolution of T μν kin
We consider first the time evolution of the energy-momentum tensor T μν kin for azimuthally symmetric initial conditions. This corresponds to solving the Boltzmann transport equation (A21) for the zeroth harmonic F 0 with initial conditions (11) brought into the form similar to (A56). The solution defines the energy momentum tensor (A13), which we shall denote with the subscript "kin" to make clear that this is the result of kinetic theory.
We start by characterizing in Fig. 2 the time evolution of the local energy density ε(t, r ) = u μ (t, r ) T μν kin (t, r ) u ν (t, r ) of a system initialized with the radial Woods-Saxon distribution (11). To explain the main features of this figure, we note the following: For a free-streaming gas of massless particles, all components propagate unperturbed along the light cone. As initially v z = 0, the particles can reach positions close to r = 0 at late times only if they originate at large radius and propagate inwards. The energy density ε(t, r = 0) at the center can therefore take non-negligible values only at times t < r tail R limited by the extent of the tail of the initial radial distribution. These particles will propagate in the (r, t) diagram on a light-like trajectory below (r, r + r tail ). On the other hand, particles that start initially at the same radial position r tail but propagate exactly outwards will show up in Fig. 2 along (r, t) = (r + r tail , r ) and all particles produced initially at smaller r or emitted not exactly inwards will stay above that line. A free-streaming density distribution is therefore a density distribution that takes non-vanishing values only in the strip limited by (r, t) = (r, r + r tail ) and (r, t) = (r + r tail , r ) in the (r, t)-diagram. This is clearly seen for the free-streaming limitγ = 0 in Fig. 2.
As one endows the system with a small interaction rate, large-angle scatterings can deflect freely propagating particles towards the neighborhood of r = 0. As a consequence, a non-vanishing energy-density around r = 0 persists for a longer time. This effect is sizeable even if scattering is infrequent: the caseγ = 0.25 in Fig. 2 corresponds to a system size that extends only over one quarter of the in-medium path length, so most particles will escape unscattered from the collision region. Yet, compared to the free-streaming case, the system maintains at its origin a non-vanishing energy density for a significantly longer time.
As one increasesγ and thus the collision propability per particle per system size, one finds that the general features described above evolve gradually, and an increasing fraction of the total energy density remains for an increasing time duration in the neighborhood of the origin of the system, see The local energy density ε(t, r ) in the (r, t)-plane, calculated from the transport equation (A54) with longitudinally boost-invariant and azimuthally symmetric initial conditions of Woods-Saxon type (11). The panels scan the only model parameter, the opacityγ , from the free-streaming limit (γ = 0) up to values ofγ sufficiently large to realize almost perfect fluid dynamic behavior. The contours depict isotherms, with ε 1/4 decreasing from one isothermal to the next by 10% of its value along the innermost contour. White arrows denote the direction u μ = (u τ , u r , 0, 0) of the radial flow field Fig. 2. In this sense, the evolution of the energy density moves with increasingγ further away from the light cone, although particles located in the tails of the system and propagating outwards will always escape along the light cone. We want to understand as a function of opacityγ , to what extent fluid dynamics provides a good description of kinetic theory. To this end, we quantify differences between kinetic theory and viscous fluid dynamics by comparing the energymomentum tensor T μν kin calculated in kinetic theory to the energy momentum tensor obtained in fluid dynamics, As we use the Landau matching condition u μ T μν kin = −ε u ν to define energy density ε and flow u μ in the kinetic theory, the ideal part T μν id = (ε + p) u μ u ν + p g μν is by construction the same, and differences between kinetic theory and fluid dynamics arise only from the shear-viscous tensor μν . In a fluid dynamic gradient expansion, μν hyd is expressed in terms of ε and u μ via the constitutive hydrody-namic relation [77] μν μν Here, the brackets around indices denote traceless, symmetrized second rank tensors, μν αβ A αβ . The terms of first and second order in fluid dynamic gradients are denoted by μν hyd 1st and μν hyd 2nd , respectively. The projector on the subspace orthogonal to the flow vector u μ is μν ≡ u μ u ν + g μν , and to first order in gradients, μν hyd is proportional to the tensor The first and second order hydrodynamic coefficients in (16) have been computed elsewhere (see, e.g. Eq. (11) in Ref. [66]) and they read Locally in space and time, we characterize the quality Q 2 (t, r ) of the agreement between fluid dynamic constitutive relation and kinetic theory by the quantity where all fields on the right hand side are understood as functions of t and r . Since . The subscript 2 indicates that Q 2 (t, r ) compares the energy-momentum tensor of the kinetic theory with the constitutive hydrodynamic equations of motion to second order in gradient expansion. It is of interest to compare also to the corresponding first order gradient expansion, and one can further consider the quantity obtained from comparing T μν kin to the zeroth order in fluid dynamic gradients, i.e., to T μν id The quantity Q 0 (t, r ) is of interest since it informs us about the extent to which the system is locally isotropic at (t, r ) in its local rest frame. The results shown for Q 0 in Fig. 3 indicate that for any non-transparent system withγ > 0, there is always a region of space-time in which the system is almost isotropic. For r = 0, the corresponding blue regions in which Q 0 is small are centered around t/R = 2. The time t/R = 2 is special in that it is the maximum time up to which particles produced initially at large radius can reach the center r = 0 along free-streaming inward-going trajectories. For smallγ , these particles will cross each other with negligible isotropizing interactions. In the local rest frame defined by Landau matching, such crossing particle streams lead to transient, locally isotropic systems. This effect is known from other fields in which one uses fluid dynamics to describe collective phenomena that is inherently non-fluid dynamical. For instance, this feature is referred to as shellcrossing in the literature on cosmological Large Scale Structure [78,79]. There, one refers to a velocity field as multivalued if the velocity of individual particles in the neighborhood of a space-time point (x, t) is not characterized by the (thermal) dispersion around the unique rest frame velocity u μ (t, x), but where it is determined instead by the (multiple) directions in which particles were located in the distant past and from which they were streaming towards (t, x). Kinetic theory realizes this scenario for almost transparent systems with small opacityγ . Non-interacting or weakly interacting streams of particles that cross each other can yield transient approximate isotropization in a neighborhood of (t, x). The tell-tale sign of this phenomenon is a negligible Q 0 (t, x) combined with a sizeable Q 1 (t, x) or Q 2 (t, x) that indicate that fluid dynamics predicts larger shear viscous components in the energy momentum tensor that is realized in the kinetic theory.
Fluid dynamics becomes a quantitatively controlled description of the collective dynamics in space-time regions where increasing the order of the fluid dynamic constitutive relation leads to a better description of the full energymomentum tensor. We therefore refer to a system evolved with kinetic equations as behaving fluid-dynamical in a space-time region around (t, r ) if Q 1 (t, r ) is sufficiently small and if this smallness persists to higher order in gradients, i.e., for Q 2 (r, t). To be specific, we use the criterion Q i (t, r ) < 0.1 which corresponds to a situation in which the summed deviations of the kinetic shear viscous tensor from its constitutive hydrodynamical form is less than 10% of the total enthalpy ε+ p of the system. This condition corresponds to the blue regions in Fig. 3. Forγ = 1, even though a small band with Q 0 (r, t) < 0.1 indicates that the system passes locally through approximately isotropic distributions, the region in which Q 1 (r, t) < 0.1 is vastly different from the region in which Q 0 (r, t) < 0.1 or Q 2 (r, t) < 0.1. This indicates that the gradient expansion does not converge and that the system displays characteristics that are very different from those of a fluid. We draw a similar conclusion for systems characterized byγ = 2, since there is no extended region in (t, r ) for which the difference Q 1 (t, r ) is small and for which this smallness persists in Q 2 (t, r ). Forγ > 4, however, our calculations indicate that the fluid dynamic gradient expansion is reliable, since the differences between T μν kin and the fluid dynamic ansatz for the energy-momentum tensor improves order by order in the  (22) gradient expansion. Further increasingγ to 8, we find that the regions of small Q 1 (r, t) and Q 2 (r, t) widen, and that they extend in particular to earlier time. In the phenomenological practice, such behavior would allow one to switch at an earlier time from a full kinetic theory description to a fluiddynamic one. For a discussion in the present framework, see, [75].
In summary, for the kinetic theory considered here, we know analytically the depth of the particle-like cut ω cut = −i/τ iso = −i γ ε 1/4 and the position of the (shear mode) pole ω pole = −i 3η 4ε k 2 in the complex plane of Fig. 1. Parametrically, an excitation of wavelength 1/k is expected to propagate more fluid-like or more particle-like, depending on whether ω cut or ω pole is closer to the real axis. Since cut-and pole-terms depend inversely on γ , it is clear that for decreasing γ , particle-like excitations dominate. Also, since we consider the evolution of the components of the energy-momentum tensor which receives contributions from all wavelengths 1/k, it is clear that the transition from a system dominated by particle-like excitations to a system dominated by fluid dynamic excitations proceeds gradually with increasing opacityγ . Quoting a precise value ofγ for this transition would rely on agreeing on a criterion (such as Q i < 0.1) of what is meant by small. Qualitatively, however, we observe from Fig. 3 that systems withγ < 2 display characteristics in the evolution of their energy-momentum tensor that show marked differences from a fluid dynamic evolution, while systems withγ > 4 evolve fluid-like from early times t/R 1 onwards. In the first case, the fluid-dynamic gradient expansion does not converge in any extended space-time region relevant for building up elliptic flow, while it does in the second case. This indicates that as the system becomes more opaque, particle-like excitations cease to dominate the collective dynamics in the range 2 γ 4.

Work in kinetic theory
The time-evolution of the transverse energy per unit rapidity and its azimuthal distribution d E ⊥ (t) dη s dφ can be calculated from kinetic theory according to (A18). The late time limit of these quantities corresponds to the experimentally accessible transverse energy distribution. In the phenomenological discussion in Sect. 5, we need to know to what extent the azimuthal average d E ⊥ (t→∞) dη s at late time differs from that at early times which is tantamount to asking how much work the system does as a function of opacityγ .
To address this question, we determine d E ⊥,free dη s by free streaming the initial conditions to late times, and by comparing to the late time behavior of an opaque system d E ⊥ (t) dη s at finiteγ . The free-streamed transverse energy d E ⊥,free dη s follows from inserting the initial condition for F(τ 0 , r, θ; φ, v z ) into (A18), and in the case of very anisotropic initial condi- tions F ∝ δ(v z ) the transverse energy remains unchanged by the free-streaming evolution. In particular, for the Gaus- This is equivalent to Bjorken's estimate of the initial energy density, Analogously, after having calculated for a system with given γ the time dependence of F(τ, r, θ; φ, v z ), we can determine dη s from (A18). In Fig. 4, we plot the resulting ratio To first order in a perturbative expansion ofγ and for Gaussian initial conditions, we know that f work (γ ) = −0.210γ (red dashed line) [73]. From the numerical results in Fig 4, we see that rare occasional collisions in almost transparent systems are most efficient in increasing the work. As collisions between particle-like excitations become more frequent with increasing opacityγ and as fluid dynamic excitations become gradually more important for the collective dynamics and the system does more work. At large opacityγ 1, the qualitativeγ -dependence of f work can be understood in a simple toy model that neglects radial dynamics: We assume that the system is free-streaming until the time τ * = τ iso (τ * , R = 0) which is the first time the system has had time to isotropize. The isotropized sys-tem then evolves hydrodynamically until it disintegrates at time τ ∼ R. The free-streaming evolution-with vanishing longitudinal pressure P L and with fixed ε(τ )τ -does not do longitudinal work and all the work is done during the hydrodynamical phase between τ * < τ < R, during which the energy density evolves with constant ε(τ )τ 4/3 . The f work is then determined by the length of the hydrodynamical phase. Solving self-consistently for the isotropization time gives which in part gives This approximate scaling for f work is observed for opacitieŝ γ 1 in the inset of Fig. 4, which shows the flattening of the rescaled quantityγ 4/9 f work (γ ) at largeγ -values.
Finally, we note the remarkable insensitivity of f work to the initial profile illustrated by the numerical similarity between the results in obtained using either gaussian (black line) or Woods-Saxon (green dashed line) initial profile in Fig. 4.

Results for the linear response v n / n
The produced transverse energy is of particular interest as it is insensitive to the modeling of hadronization. It is defined by the first p ⊥ -moment of the measured particle distributions The Fourier coefficients v n and reaction plane orientations ψ n parametrize the azimuthal dependence of the transverse energy flow. Here we have used the fact that particles with momentum rapidity y after the last interaction will end up in a spacetime rapidity η s = y. See Appendix A2 for details. Spatial deformations δ n,m are introduced in the initial conditions according to (8). For the kinetic theory formulation explored in this work, the first results for the linear response of the elliptic energy flow v 2 to initial elliptic perturbations 2 were reported in Ref. [75]. The purpose of the present subsection is to complete this information with results on higher harmonics, and to demonstrate that the results depend negligibly on the transverse profile function. These results are complete to all orders inγ , but they are first order in n . The formulation in Sect. 3 one allows for an extension to all orders in n in future work.
This standard definition of n is adopted for all initial conditions irrespective of their radial dependence (∝ δ n,m r m ) in (8) and their azimuthally averaged transverse profile. The black and the purple lines in Fig. 5 correspond to Gaussian background initial conditions of (10) and (8) with δ 2,2 = 0 (black line) or δ 2,3 = 0 (purple line). This illustrates the insensitivity of the response to details of the radial profile of the perturbation. Moreover, there is a remarkable insensitivity to the background profile, since Wood-Saxon (yellow line) and Gaussian (black line) initial conditions yield numerically similar linear responses. Furthermore, the green line in Fig. 5 corresponds to an isotropic initial momentum distribution F(r, θ, v z ) = ε 0 e −r 2 /R 2 with a δ 2,2 perturbation. Again no significant sensitivity to the initial condition is visible. Similar observations can be made for v 3 / 3 from the bottom panel of Fig. 5. Fig. 5 illustrates for an important class of observables how kinetic theory interpolates between the limits of freestreaming and perfect fluidity. In the free streaming limit (γ = 0) momentum asymmetries v n are not built up from initial spatial asymmetries n . As a consequence, v n / n starts from zero, and itsγ -slope (red lines in Fig. 5) close toγ = 0 can be understood as the consequence of perturbatively rare final state interactions in a small and/or dilute system. In the opposite limit of ideal fluid dynamics (γ → ∞), momentum asymmetries are built up from spatial eccentricities with vanishing dissipation and thus most efficiently. As a consequence, v n / n increases monotonously with increasing opacityγ and it approaches at very largeγ the ideal fluid limit (dashed blue lines in Fig. 5) from below. Forγ < 16, the difference between this ideal fluid limit and the result of full kinetic theory indicate dissipative effects. As discussed in the context of Fig. 3, these dissipative effects are accounted for by viscous fluid dynamics or, for smallerγ , by particle-like excitations outside the fluid dynamic description.

Data comparison
The kinetic theory defined in Sect. 3 one and analyzed in Sect. 4 is a one-parameter model. By comparing it to experimental data, we do not intend to claim that a complete phenomenological understanding of data could be based on this one-parameter model. It is obvious that it cannot, and by stating at the beginning of this section the obvious, we hope to preempt any possible misunderstanding of this point. Rather, by comparing to data, we want to understand to what extent the one generic physics effect implemented in this model could account for one central qualitative feature in the data, namely for the remarkable system size dependence of measures of collectivity. Is the centrality dependence of flow measurements consistent with a dynamics that interpolates between free-streaming and viscous fluid dynamics, or is it more consistent with a collective dynamics that is always fluid-like in the sense that it depends negligibly on particlelike excitations even for the smallest collision systems studied at the LHC?

The centrality dependence of opacityγ
We recall that in kinetic theory, the retarded two-point correlation functions of the energy-momentum tensor display fluid-dynamic poles and particle-like cuts. The position of the hydrodyamic poles is unambigously determined by the parameters of the kinetic theory. In particular, for the kinetic follows from Eq. (18) if one choses ε = 13 T 4 consistent with results from lattice QCD calculations [80,81]. As discussed in Sect. 3 one, the numerical solutions of this kinetic theory depend only on the dimensionless opacitŷ γ = γ R 3/4 (ε 0 τ 0 ) 1/4 , which we rewrite with the help of Eqs. (24), (26) and (31) in the form Here, the factor f work (γ ) appears since we have traded the initial energy density ε 0 for an expression in terms of the experimentally measured final transverse energy d E ⊥ dη s . The function f work (γ ) is a prediction of kinetic theory.
By solving (32) forγ , we can determine the centralitydependence ofγ from the centrality dependence of the transverse energy d E ⊥ dη s and the centrality dependence of the rms radius R of the initial transverse energy profile. We note that for LHC energies, the transverse energy d E ⊥ /dη varies by a factor 80 between the most central (1%) and most peripheral (80%) centrality, while the rms radius R changes by less than a factor 2. As a consequence, the centrality dependence ofγ is largely determined by the centrality dependence of the transverse energy, and modeling uncertainties of R play a minor effect in determiningγ (e.g., varying R by 10% amounts only to doubling the line width in Fig. 6). This allows us to determine for simplicity the rms radius R from an optical Glauber model [82].

Opacityγ in PbPb at the LHC
The centrality dependence ofγ shown in Fig. 6   values determined from the p ⊥ -spectrum d N dp ⊥ dη s published by ALICE in Ref. [83] for nine centrality bins between 0-5% and 70-80%, corrected for the presence of neutrals and interpolated to finer intermediate centrality bins. Combined with the conclusions drawn from Fig. 3, Fig. 6 informs us for which values of η/s fluid dynamic behavior can be expected as a function of system size.
For instance, for η/s = 2/4π , the valueγ ≈ 8 reached in the most central PbPb collisions (see, Fig. 6) signals according to Fig 3 that for matter of this viscosity the collective dynamics realized in central PbPb collisions is fluid-like. But for semi-peripheral PbPb collisions with centrality > 50%, the η/s = 2-curve in Fig. 6 shows valuesγ < 4 for which a fluid dynamic picture of the evolution becomes questionable, and for > 80% centrality, valuesγ < 2 signal a collective dynamics in peripheral collisions which is qualitatively different from that of a fluid. As seen from Fig. 6, the centrality range for which a fluid dynamic gradient expansion is expected to be quantitatively reliable (γ > 4) or qualitatively indicative (γ > 2) changes rapidly with η/s. According to Fig. 6, a gradual numerical change form η/s = 1/4π to η/s = 4/4π translates into a qualitatively different picture of the physics at work, since it translates into scenarios in which particle-like excitations make either the negligible or the dominant contribution to the collective dynamics in semi-peripheral and peripheral PbPb collisions at the LHC. Irrespective of the specific kinetic theory studied in the present work, this provides a strong physics argument for a precise determination of η/s from data.

Opacityγ in AuAu at RHIC
In close analogy to the calculation ofγ for PbPb collisions at the LHC, we calculateγ for √ s NN = 200 GeV AuAu collisions based on STAR data of the N part -dependence of d E ⊥ /dη s [84]. For the same centrality and the same assumed value of η/s,γ is approximately 20% smaller at RHIC than at the LHC, and the parameter regions in which fluid-like and particle-like excitations dominate collectivity shift accordingly, see, Figs. 6 and 7. For nucleus-nucleus collisions at RHIC to show a more fluid-like behavior (i.e. correspond to a largerγ ) than at the LHC, the matter produced at RHIC would thus have to display a significantly smaller η/s, since at comparable values of η/s, the higher center of mass energy and the higher transverse energy density attained at LHC facilitates fluid-like behavior.

Opacityγ in pPb at the LHC
To calculate the opacityγ for pPb collisions at the LHC, we make the simplifying assymption that the rms-width R of the transverse area initially active in pPb collisions does not change with event multiplicity. This assumption has a physical and a pragmatic motivation. Physically, it amounts to stating that the dispersion in the multiplicity distribution at fixed impact parameter is so broad that the correlation between event multiplicity and event centrality (i.e., impact parameter) is weak. From a pragmatic point, choosing R to be multiplicity-independent is reasonable, since we do not have firm knowledge about how the initial geometry of pPb collisions changes with multiplicity. Moreover, for the purpose of calculatingγ , this lack of knowledge is partially alleviated by the fact that information about the multiplicity-dependence of R entersγ only as a fourth root. We therefore choose in the following R = 1 fm as default value for pPb collisions, noting that even extreme values of R = 2 fm would affect results forγ only by 20%.
With the assumption of a multiplicity-independent geometry in pPb, the entire centrality dependence ofγ in Eq. (32) originates now from the experimentally inferred N partdependence of d E ⊥ /dη s that we take from measurements of the CMS collaboration [85]. The results of Fig. 8 show a substantial reduction ofγ in the smaller pPb collision systems if compared to heavy ion collisions. In pPb collisions, a transport coefficient η/s = 1/4π indicating maximal fluidity needs to be assumed to reach at least in the highest multiplicity-bins values ofγ > 4. In summary, under essentially all conditions explored in Fig. 8, value ofγ is so small in pPb collisions that collective dynamics cannot develop fluid-like behavior.

Comparing kinetic theory to v n / n (γ ): assumptions and uncertainties
Before presenting in the following subsections comparisons of kinetic theory to data on v n / n , we list here the main elements of this comparison and we discuss the underlying assumptions and resulting uncertainties. While data on d E ⊥ dη s are published for all collision systems, data on the harmonic decomposion (29) of d E ⊥ dη s dφ in terms of energy flow coefficients v 2 and v 3 are measurable but have not been published yet. Therefore, we have to infer them from the measured particle flows v particle m ( p ⊥ ), m = 2, 3 and from the published single inclusive charged hadron spectra where d E ⊥ dη s = 3 2 ∞ 0 dp ⊥ p ⊥ d N dp ⊥ dη s (the factor 3 2 corrects approximately for neutrals). The integrands in (33) differ by one power of p ⊥ from the ones used to determine the p ⊥integrated particle flow Data on v particle m are usually quoted as p ⊥ -integrated within a finite p ⊥ -range that varies between experiments. In all collision systems discussed in the following, we have determined v energy 2,3 from v particle 2,3 and d N dp ⊥ dη s by first checking consistency of the published p ⊥ -integrated particle flows with (34), and then determining the energy flow coefficients according to (33). One can then determine the correction factor c e−corr between measured particle flow and extracted energy flow, Both particle flow and momentum flow coefficients (33) and (34), respectively, are centrality dependent and therefore, in principle, c e−corr,m could be centrality-dependent as well. We anticipate that the centrality-dependence of c e−corr,m turns out to be negligible in all collision systems considered. Particle flow coefficients v particle m are measured with different methods that are designed to include different sources of particle correlations and that therefore yield numerically different results. In particular, the 2nd and higher order cumulant flows v particle m .., are designed such that they include only particle correlations that persist in connected 2-, 4-, 6-, ... particle correlation functions. Moreover, to eliminate contributions from jet-like structures, data for v particle m are often designed to include only the correlation of particles separated by a sizeable rapidity gap | η s |.
To choose the variant of particle flow measurements from which we infer v energy m via (33), we start from the following two conisiderations. First, one of the simplest perturbative corrections to free-streaming that is accounted for by kinetic theory may be thought of as a single rare scattering that leads to correlations amongst very few of the final state particles and that leads, a fortiori, to an azimuthal correlation of d E ⊥ dη s dφ carried by a small fraction of the entire d E ⊥ dη s . As such a correlation is not shared by all particles in the event, it would not fully persist in higher order cumulants, while it is fully accounted for in kinetic theory. Second, the formulation of the kinetic theory presented here is not sufficiently detailed to include the characteristic dynamics of jet-like fragmentation patterns, and a meaningful comparison of kinetic theory to data should thus start from an energy flow v energy m from which jet-like correlations are eliminated to the extent possible. The combination of both considerations prompts us to calculate v energy m from 2-particle cumulant measurements with rapidity gap. Consistent with this choice, we use throughout this work second order cumulants to characterize spatial eccentricities, Based on these considerations, the following data sets enter our analyses: 1. Data for PbPb collisions at the LHC p ⊥ -differential particle flow v 2,3 {2, | η s |}( p ⊥ ) [2] published in nine centrality classes between 0-5% and 70-80%, and single inclusive charged hadron spectra [83] in the same centrality classes and the same mid-rapidity range are used to determine c e−corr,2 = 1.34(1.33) and c e−corr,3 = 1.52(1.51) for √ s N N = 5.02 ( √ s N N = 2.76) TeV PbPb collisions. Fluctuations of these numbers with centrality are small, do not reveal any systematic trend, and may be attributed to uncertainties of p ⊥differential flow data at higher p ⊥ . The p ⊥ -integrated particle flow calculated from these data according to (34) is checked to be consistent with the p ⊥ -integrated v 2 {2, | η s | > 1} [2] restricted to the experimentally chosen range 0.2GeV < p ⊥ < 3GeV. The factors c e−corr,2 and c e−corr,3 are then used to infer the energy flow via (35) for the 80 different 1%-wide centrality bins published in [2] for √ s N N = 5.02 TeV PbPb, and for the somewhat fewer bins between 0% and 80% centrality published for √ s N N = 2.76 TeV PbPb. The data on single inclusive charged hadron spectra [83] are further used to determine d E ⊥ /dη s for the nine bins between 0% and 80% centrality, and are then interpolated to obtain centrality dependent values for d E ⊥ /dη s in 80 1%-wide centrality bins. This is checked to be consistent with published results on the centrality dependence of d E ⊥ /dη s . In addition to these data of the ALICE collaboration, there are CMS-data on v 2 {2, | η s | > 2} in √ s N N = 5.02 PbPb [4] that could allow one to extend particle flow measurements to ≈ 93% centrality (Ref. [4] reports data on pp and pPb, but its HEPDATA-repository documents the latest CMS data on PbPb collisions). As these data are presented as a function of off-line tracks N offline trk , as they are integrated over a slightly different p ⊥ -range, and as they are taken in a much wider rapidity range, there are subtelties in consistently infering energy flow as a function of centrality with CMS and ALICE-data on the same plot. We are in contact with both collaborations to arrive at such a combination, but this lies outside the scope of the present manuscript.

Data on pPb collisions at the LHC
CMS data on the pseudorapidity dependence of the transverse energy d E ⊥ /dη s in pPb collisions at √ s NN = 5.02 were presented in [85] as a function of N part , and pPb particle flow measurements v 2 {2, | η s | > 2} were presented in [4] against the experiment-specific number of offline tracks N offline trk . Information tabulated in Ref. [86] relates N offline trk to the fraction of the total pPb cross section (i.e. centrality), and this average centrality can be converted into an average N part , as done for instance in Ref. [87]. However, CMS data points for N track offline 150 that correspond to the < 1% most active events of the total cross section cannot be related meaningfully to N part .
As v 2 {2, | η s | > 2} for N track offline 150 is approximately constant [85] information from all data points at larger N offline trk can be regarded as grouped into the data point of largest N part . In this way, both particle flow and transverse energy in pPb at the LHC is known as a function of N part . The following discussion uses the N part -dependence of v 2 {2, | η s | > 2} and d E ⊥ /dη s thus obtained without relating N part to a geometrical picture of the collision. The same factor c e−corr,2 = 1.34 determined for √ s NN = 5.02 PbPb collisions is used to determine energy flow, since the absence of sufficiently fine-binned p ⊥ -differential v 2 data does not allow one to perform the same determination of c e−corr,2 = 1.34 as for PbPb. We anticipate that none of our conclusions will depend on the precise numerical value of c e−corr,2 = 1.34 used in pPb.

Data for AuAu collisions at RHIC
The centrality dependence of d E ⊥ /dη s is taken from Ref. [84]. We infer energy flow v 2 for √ s NN = 200 GeV Au+Au collisions at RHIC by paralleling closely the procedure described for PbPb collisions at the LHC. STAR p ⊥ -differential charged hadron v 2 ( p ⊥ ) data [88] in the wide 10-40% centrality bin is weighted with charged single inclusive hadron spectra from [89] according to Eqs. (33) and (34) to infer a correction factor c e−corr,2 = 1.5 for p min ⊥ = 0.1 GeV. This factor is applied to particle flow data of p ⊥ -integrated v 2 (| η s |) presented in 11 N part -bins by PHOBOS [7]. To allow for a better comparison with other calculations presented in this paper, the N part -dependence is then converted into the impact parameter dependence obtained from an optical Glauber calculation [82].

Non-ideal equation of state
The kinetic transport theory studied here is for a conformal system with c 2 s = 1/3. While this assumption may be relaxed in future work, we have not calculated yet v n / n (γ ) for more realistic equations of state in kinetic theory. Recent lattice QCD results [80,81] indicate values c 2 s (T ) which decrease gradually with decreasing temperature, reaching values 0.23, 0.25 and 0.29 at temperatures T ≈ 200, 225 and 300 MeV, respectively. A more realistic formulation of the expansion would dynamically average over this temperaturedependence. To estimate the difference that such a formulation could make, we calculate first from ideal fluid dynamics the correction factor c eos with which the ideal hydro baselines (blue dashed lines in Fig. 5) must be multiplied to reflect a reduced sound velocity c 2 s . We find that in ideal fluid dynamics with c eos = 0.86 (0.93) for c 2 s = 0.25 (=0.29), respectively. For sufficiently largeγ , the kinetic theory result for v n / n approaches this ideal hydro baseline from below, and must be affected by the same factor c eos in case of a non-ideal equation of state. To estimate uncertainties, we now assume that this correction factor c eos can be applied to all values ofγ . We choose in the following c eos = 0.9 as a default (corresponding to c 2 s = 0.27), and we compare to c eos = 0.8, which-according to these numerical estimates-may be viewed as a conservative overestimation of the expected correction. Fig. 9 The centrality dependence of the (second order cumulant) elliptic eccentricity 2 = 2 {2} (upper plot) and triangular eccentricity (lower plot) for two models discussed in the text. Data taken from Ref. [92] and displayed over a wider range of impact parameter. The grey band is used to characterize the model-uncertainty of these eccentricities in the following study

Non-linear response to spatial eccentricities
The kinetic theory calculations of v n / n (γ ) in Fig. 5 determine only the linear response v n to an initial spatial eccentricity n Both theoretical and experimental evidence indicates that the true response has a subleading non-linear component that is non-negligible in particular for v 2 in mid-central collisions where eccentricities are large. Future calculations may allow us to determine these non-linear corrections. Here we estimate them by assuming that they have the same size in kinetic theory as in fluid dynamic simulations in which they have been quantified. To this end, we recall first that symmetry arguments dictate that the leading non-linear corrections to v 2 is O 3 2 . From fitting the visible nonlinearities in the Fig. 17 of the fluid dynamic study [90], we find v particle 2 Non-linear corrections have been determined also in Fig. 2 of Ref. [91]. According to this study, c nl ≈ 0.6-0.65 for LHC PbPb collisions at 50% centrality while more central collisions show a significantly smaller c nl (c nl ≈ 0.3-0.35 for 5% centrality). However, the precise value of c nl in central collisions is unimportant for our purpose since 2 is small and non-linear effects are known to be negligible. To estimate uncertainties of the non-linear response c nl , we therefore compare in the following the case of a negligible nonlinearity (c nl = 0) with the case of the maximal non-linearity (c nl = 0.75) supported by a study. The same correction is also applied to v 3 , where non-linear corrections are much less important since 3 is smaller.

Centrality dependence of initial eccentricities 2 and
3 In comparing data to the kinetic theory prediction of v n / n (γ ), the model-dependence of the initial eccentricities n needs to be accounted for. Here, we consider two recent initial state models that were studied in Ref. [92] (referred to as GGMLO in the following) and that are based on different physics assumptions. They compare one of the current phenomenological standards, the so-called Trento model for p = 0 [93], to a model based on saturation physics for which the variation of the saturation scale Q s translates into a band in the centrality dependence of initial eccentricities. Ref. [92] shows the corresponding eccentricities for impact parameter b < 10 fm, since saturation physics is thought to lose predictive power for more peripheral collisions. As we study in the following systems up to 80% centrality, we require the impact parameter dependence of n for b 13.9 fm. We thank the authors of Ref. [92] for providing their data for this extended regime which we replot in Fig. 9. We have further obtained corresponding data sets for 2.76 TeV PbPb collisions that take for the saturation model slightly larger values of the eccentricity (data not shown) and that we shall use. The gray band in these figures will enter as model-dependent uncertainty in the following study.

Confronting kinetic transport to data on v n / n (γ ) in
PbPb at the LHC We can determine now for each centrality bin in each collision system i. the energy flow v n from data, ii. the transverse energy per unit rapidity d E ⊥ dη s from data, iii. the rms radius R of the initial distribution as a function of centrality from a Glauber model, iv. the eccentricity n .
How to arrive at these quantities has been detailed in the previous subsection. The resulting basis for comparing data to the kinetic transport formulated in Sect. 3 one can be summarized in the simple expression Here, the quantities c e−corr , c eos and f corr are fixed as discussed in Sects. 5.2.1, 5.2.2 and 5.2.3, respectively, and the Fig. 11 Same as Fig. 10, but for the ratio v 3 / 3 as a function ofγ . Since 3 is significantly smaller than 2 , correcting or not correcting with Eq. (39) for the non-linear 3 -dependence of v 3 makes a negligible effect model-dependent eccentricities 2 , 3 are varied according to Fig. 9. The first line (39) represents the centrality-dependent experimental input. The factor ∂v n ∂ n in the second line (40) represents the theoretical expectation. It is calculated from kinetic transport and its centrality dependence arises solely from itsγ -dependence.
We use data for 80 1%-wide centrality bins in 5.02 TeV PbPb collisions to evaluate (39) and to plot the results against the 80γ -values of the corresponding centrality bins (see Eq. (32) and Fig. 6). Figures. 10 and 11 show the comparison of these data to kinetic theory results for v 2 2 corr (γ ) and v 3 3 corr

(γ ).
As seen from Fig. 6, the 80 bins between 1% and 80% centrality correspond to ranges 4 γ 16 ( 2 γ 7.5, 1.3 γ 5.0) for η/s = 1/4π (η/s = 2/4π , η/s = 3/4π ), respectively. Accordingly, the blue, red and brown bands for η/s = 1/4π (η/s = 2/4π , η/s = 3/4π ) in Fig. 10 scan over the correspondingγ -ranges, with the most peripheral centrality bin determining the left-most data of the corresponding band, and the most central bin determining the right-most data point. For fixed value ofγ , the sizes of the plotted error bars reflect only the modelling uncertainties in the initial eccentricities which are depicted by the grey band in Fig. 9. These uncertainties get larger for more peripheral centralities, when the uncertainties in the models (see, Fig. 9) are larger. The uncertainties are also somewhat larger in the most central bins, when v 2 2 corr is obtained by dividing v 2 with a small eccentricity 2 of non-negligible uncertainty.
To illustrate the size of the correction terms in Eqs. (39), (40), we display the same data in the lower panel of Fig. 10 without correcting for nonlinearities, that is, by setting f corr = 1. Similarly we expose the theoretical uncertainty arising from the equation-of-state correction factor c eos by We find that the centrality dependence of v 2 2 corr extracted from data coincides well with that from kinetic theory for η/s ≈ 3/4π if the corrections c eos = 0.9 for the equation of state and c nl = 0.75 for the non-linear response are applied. For the conservative overestimate c eos = 0.8, or for a significantly smaller non-linear correction c nl , a value close to η/s = 2/4π would be favored. In contrast, neglecting a correction for the non-conformal equation of state (c eos = 1.) would result in a favored value η s ≈ 4/4π . A more precise knowledge of c eos and c nl may allow one to provide a tighter constraint on η/s, but within the above-mentioned sizeable uncertainties, we restrict ourselves to the conclusion that the favored value of η/s is located in the range 2/4π < η/s < 4/4π . To avoid duplicating all subsequent figures, we show and discuss from now on only the case c nl = 0.75. For triangular flow v 3 , the non-linear correction is comparatively small since 3 is smaller, and for elliptic flow, the effect of the correction c nl = 0.75 is comparable to that seen in Fig. 10.
This analysis is repeated for triangular flow in Fig. 11. Since the maximal 3 values are significantly smaller than  Fig. 9. Correcting for a QCD-like equation of state with c eos = 0.9, a value η/s = 3/4π is favored those for 2 (see, Fig. 9), a correction for a non-linear dependence on eccentricity f corr ( 3 ) turns out to be negligible (data not shown). Also, since the gray uncertainty band shown for the centrality dependence of eccentricity in Fig. 9 is in semicentral collisions much broader for 3 than for 2 , the corresponding uncertainties in plotting v 3 3 corr in Fig. 11 are more pronounced for data from semi-central events. We see that kinetic theory accounts for the magnitude and centralitydependence of v 2 2 corr and v 3 3 corr with values of η/s that lie within the same ball-park, although the triangular flow data prefer slightly lower η/s or-alternatively-prefer a slightly smaller 3 than the one displayed in Fig. 9.
The analysis of the elliptic and triangular flow are repeated for √ s NN = 2.76 TeV PbPb collisions in Fig. 12. Compared to the same centrality bin at the higher center of mass energy of 5.02 TeV, the smaller produced transverse energy results for the same centrality in slightly smaller values of γ . The conclusion about the preferred parameter range for η/s remains approximately the same for √ s NN = 2.76 and √ s NN = 5.02 TeV PbPb collisions. Given the simplicity of the kinetic theory studied here and the uncertainties enumerated in Sect. 5.2, we do not draw too tight numerical conclusions about η/s from the present study. We solely note that the parameter range 2/4π < η/s < 4/4π is favored by this analysis. We further note that this parameter range is consistent with the values favored in recent kinetic theory studies (0.2 < η/s < 0.4, [31]), and it is consistent with a full phenomenological analysis that includes a hydrodynamic evolution stage and that favors η/s = 2.5/(4π) [94]. Other full phenomenological analyses [95] favor smaller values for η/s. The constraint η/s < 4/4π is significantly lower than the value favored in our earlier analysis [73]. We note that this earlier analysis was based on the single hit line in Fig. 5, and that the difference compared to the present study arises from a combination of several improvements. In particular, the present analysis includes non-perturbative control over thê γ -dependence beyond the linear approximation, it takes into account the correction factors in (39), (40), and it exploits the full experimentally available centrality dependence.
In combination with the η/s-dependence ofγ , the loose constraint 2/4π < η/s < 4/4π implies an interesting qualitative statement. Irrespective of where the favored value of η/s is located in the range 2/4π < η/s < 4/4π , semiperipheral or peripheral PbPb collisions reach down to valueŝ γ < 2 for which a fluid dynamic picture does not correspond to the physical reality, while the most central PbPb collisions reach up to valuesγ > 4 for which a fluid dynamic picture is a suitable interpretation of the physical reality.
To further emphasize this point, Fig. 13 compares kinetic theory with η/s = 3/4π to elliptic and triangular energy flow data without deviding them by eccentricity. Centrality ranges of PbPb collision correspond-provided that the chosen values of η/s is physically realized-to ranges ofγ for which a physical picture of the collision based on fluid dynamics is (γ > 4) or is not (γ < 2) suitable. The uncertainties in favoring η/s = 3/4π can be judged from the parameter scans in Figs. 10 and 11. We note as an aside that Figs. 10 and 11 for η/s = 3/4π and Fig. 13 contain identical information; any visual impression that the agreement between theory and experiment is better in Fig. 13 is erroneous.

The centrality dependence of eccentricities 2 and 3 preferred by kinetic theory
Rather than testing a combination of kinetic theory and model-dependent initial eccentricity versus flow data, as done in the previous subsection, we invert now the logic of the analysis. We assume that the kinetic theory reflects physical reality and we use it to extract the eccentricity of the collision system from data. To this end, we rewrite (39), and we solve for n for each centrality bin. As this procedure can be followed without making assumptions about n , it allows us in the following to compare different collision systems (PbPb and pPb at the LHC, AuAu at RHIC) without Fig. 14 The eccentricities 2 and 3 preferred by kinetic theory with c eos according to (41) and compared to two eccentricity models of Ref. [92] discussing the relative uncertainties with which the initial eccentricities of these systems are thought to be known. The centrality dependence of n thus determined depends on η/s, and the physical value of γ = 0.11/(η/s) can be constrained from it to the extent to which information about the size and centrality dependence of n is assumed to be known. We start this discussion with the 5.02 TeV PbPb data, for which Fig.14 shows the centrality dependence of eccentricities 2 , 3 reconstructed from (41) and overlaid with the eccentricity models depicted in Fig. 9. The conclusion is obviously the same as the one reached in Sect. 5.3: combining LHC data and the eccentricity model of Fig. 9 favors a value of η/s close to 3/4π . In addition, (41) informs us by how much one would have to distort the model of initial eccentricity to favor a value outside the range 2/4π < η/s < 4/4π .
We have reconstructed the eccentricity 2 for 200 GeV AuAu collisions from data sets discussed in Sect. 5.2. The initial geometrical distribution in AuAu collisions at RHIC and in PbPb collisions at the LHC can differ not only because of the slighly different nuclear size, but also because particle production processes and the ensuing event-wise fluctuations are expected to evolve over an order of magnitude in √ s NN . Because of the uncertainties discussed in Sect. 5.2, the reconstructed eccentricities displayed in Fig. 15 do not exclude small difference in the eccentricity of both collision systems. Such differences could arise within the present analysis for instance if the value of η/s varies significantly with √ s NN , or if the factor c eos in (41) takes a slighly different value for both collision systems. Also, the slightly different size of the Au and Pb nucleus implies that larger differences would certainly become visible if RHIC data were extended to more peripheral collisions. Despite these caveats, Fig. 15 seems to supports an approximately √ s NN -independent elliptic eccentricity.

Reconstructing eccentricity for pPb collisions
The eccentricity in pPb collisions is particularly difficult to infer from multiplicity distributions close to mid-rapidity, since the correlation between event multiplicity and transverse geometry weakens when the dispersion of the multiplicity produced at fixed impact parameter approaches the width of the entire multiplicity distribution. In plottingγ for pPb collisions in Fig. 8, we had assumed already the limiting case that this correlation is negligible so that events of different multiplicity (classified experimentally as a function of N part ) correspond to systems of initially similar transverse width R.
In reconstructing eccentricity according to (41) from the 5.02 TeV pPb data discussed in Sect. 5.2, Fig. 16 gives independent support to this assumption of a multiplicityindependent initial geometry of the active area in pPb col- The flat eccentricity in Fig. 16 is qualitatively different from the one predicted by Glauber-type models for which pPb 2 decreases with increasing N part [96]. These Glaubertype models have been used to initialize phenomenologically successful fluid dynamic simulations of pPb collisions [96,97]. In this way, these studies support the view that pPb collisions can be viewed as creating initially a small fluid droplet, in which fluid-dynamic expansion translates an impact parameter-dependent spatial anisotropy into a centrality-dependent momentum anisotropy. However, if we try to translate the geometrical picture of pPb collisions documented in [96,97] into the parameters of the present study, we systematically arrive at values ofγ that are too small to support a fluid dynamic picture. Rather, insensitive to initial geometry, the kinetic theory studied here builds up flow under spatio-temporal constraints in which fluid-dynamic modes cannot dominate collectivity, and the multiplicitydependence of v 2 / 2 can be viewed as arising solely from the multiplicity-dependent scattering probability that enhances collectivity within a system of fixed eccentricity.

Comment on measures of system size in small systems
As explained in Sect. 3 one, the opacityγ is the only dimensionless variable that observables may depend on, and hence it is the unique variable in which the system size can be measured. It is natural to think of system size in units of mean free path. In the present model, this quantity at the relevant time when flow is built up (τ = R) at r = 0 reads This quantity can be calculated exactly in the current kinetic theory setup. Denoting ε(τ =R)R The fraction f work (γ ) is closely related to f 0→R , and in thê γ -range studied here, we find e −1 f work (γ ) = f 0→R (γ ) up to 15% corrections. Therefore, R τ R ≈γ up to small corrections in large systems, for which the scaling of Eq. (28) suggests that this quantity is proportional to Like the opacityγ used in the plots and defined in Eq. (32), the ratio R τ R (τ =R) is a suitable scaling variable in the sense that it depends only onγ and that it increases monotonically withγ .
The energy density at the origin at time τ = R may be estimated also in other ways. For instance, following Bjorken, one may assume that the system decouples at time τ = R and that it evolves free streaming from this time onwards. A simple geometrical picture yields then where ε ⊥ (R) is the transverse energy density constructed from the total transverse energy d E ⊥ /dη s within an appropriately chosen transverse area A ⊥ ∼ π R 2 . Using this estimate to evaluate τ R = 1 γ ε 1/4 , one finds where a = ε ⊥ /ε = [ 2 π , 1]. Here, the value a = 1 is realized for a maximally anisotropic system ( p z = 0) at time τ = R and a = 2 π for an isotropic one. The parametric estimate (46) is consistent with the full calculation (43).
In general, the scaling variable (43) is proportional to the fourth root of the transverse energy. It is noteworthy that for the specific case of large systems (γ 1) which had time to reach local thermal equilibrium, scaling in opacityγ becomes interchangeable with a scaling in d N/dη s . To see this, let us specify the assumption that the system has reached thermal equilibirium such that p ⊥ = f ε ε 1/4 , where f ε is a constant for a conformal system. Inserting this into the expression forγ in terms of d N dη s , and solving self-consistently forγ , one findŝ where . Equation (47) holds for systems of any size, but it is only in the limit of large opacity that the scaling variableγ seizes to depend on the geometrical extent and depends solely on event multiplicity [64]. This is consistent with the wide-spread experimental practice of using d N dη s as a proxy of system size and as a scaling variable for comparing different collision systems. So, according to this kinetic theory, for small systems, measurements will take system-independent values if plotted againstγ , while they take in general system-dependent values if plotted against d N dη s .
We finally comment on differences between our formulation of kinetic theory and formulations of massless kinetic theory with fixed cross section σ . In these latter models, the mean free path τ R (τ ) = 1 n(τ ) σ at time τ may be expressed in terms of the particle density n(τ ) ∼ 1 τ A ⊥ d N dη s . In a small, low-density system, the first correction to free-streaming is expected to be proportional to the inverse Knudsen number that may be written as , [for a small conformal system].
The scaling of Eq. (51) is that of a system with conformal symmetry, while (50) arises from a model of explicitly broken conformal symmetry. Up to corrections due to renormalization group running, conformal symmetry is realized in high energy QCD in small systems, as well as in high-temperature QCD. There are additional scales that arise at lower energies, including e.g., particle masses. The introduction of a fixed cross section may be viewed as a model of non-conformal matter.
It is important to understand whether conformal scaling (51) of anisotropic flow with multiplicity can be established experimentally in the smallest systems. Whether or not corrections to this scaling could be established, would inform us about important elements of the microscopic dynamics underlying collectivity.

Conclusions
To ask what is the microscopic structure of quark-gluon plasma, is to ask how the plasma behaves away from the hydrodynamic limit. The study of the onset of signs of collectivity as a function of system size offers one setup where this question can be addressed. The effective description of the plasma in terms of fluid dynamics must eventually break, giving way to non-hydrodynamic modes to start dominating the dynamics, and the precise way how this happens reflects the microscopic structure of the plasma. Here, we have studied the consequences of assuming that the plasma has an underlying quasiparticle description. We have chosen this as our starting point because it is consistent with free-streaming particles in small systems.
Many microscopic models may appear to be consistent with the data on small systems and given the large unknowns in the initial condition of the collision it may be difficult to discriminate amongst them. One qualitative feature that quasiparticle models have in common is that the location of the hydrodynamic pole (the value of η/s) determines also the location of the non-hydrodynamic sector in the complex frequency plane. While the details of this correspondence may vary between different specific quasiparticle models, the scale at which non-hydrodynamic structures become dominant R is given by the mean free path R ∼ 1/(γ ε 1/4 ) which also determines the specific shear viscosity η/s ∼ 1/γ . Therefore the knowledge of transport properties of the plasma also uniquely determines how the free-streaming behaviour is reached in systems that have transverse sizes smaller than the mean free path. In Sect. 5 we have seen that this prediction is not in contradiction with the available data.
Based on this exercise, to what extent can we then conclude that the plasma does have a quasiparticle description? In quantum field theory at weak coupling, there are additional structures in the complex plane to the hydrodynamic poles and the quasiparticle cut. These appear at the scale of the first Matsubara mode Im [ω] ≈ −4π T . They reflect the quantum mechanical nature and the de Broglie wavelength of the quasiparticles. When the mean free path becomes ∼ 1/T ∼ 1/ε 1/4 , the quantum field theory becomes strongly coupled and there is no clear separation between the quasiparticle cut and the quantum mechanical Matsubara modes. We have no firm knowledge about what happens in quantum field theory at these couplings, but we do know that in the limit of infinite coupling the nonhydrodynamics structures still lie at the scale of Im [ω] ∼ −T . Therefore one could expect that in a model that goes beyond quasi-particles, the scale at which non-hydrodynamical modes appear saturates such that R ∼ min( 1 γ ε 1/4 , 1 ε 1/4 ). This is qualitatively in contrast to the quasiparticle models for which R ∼ 1 γ ε 1/4 . It may be that such models also describe the data well, but given the good agreement with the data to the quasiparticle model, we must acknowledge that we do not have evidence that the quark-gluon plasma does not have quasiparticle excitations.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Author' comment: This manuscript has no associated data. All results shown in the figures follow directly and unambiguously from the formulas displayed in the manuscript.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Free-streaming variables
and which are consistent with Eqs. (A21), (A22). Since the coefficient functionsF (n) (t, x ⊥ ; φ,ṽ z ) are invariant under azimuthal orientations, they can depend on azimuthal angles only viaφ r = φ −θ . Transforming the integrand of (A25) to free-streaming coordinates with the help of Eqs. (A19) and (A16), we write We shall finally determine the linear response coefficients v n δ n by comparing the parametrization (29) to the late time limit of (A28). To this end, we solve the kinetic evolution equations (A26), (A27) forF (n) in free-streaming coordinates.

b. Isotropization time approximation (ITA)
The isotropization time approximation is based solely on the assumption that f evolves towards a distribution f iso ( p μ u μ ) which as a consequence of being isotropic can depend only on the scalar p μ u μ . In general, this assumption is not sufficient to specify the functional form of f iso . Remarkably, however, this assumption fully specifies the functional form of the first momentum moment Here, the dependence on 1/ −v μ u μ 4 follows from the dimensionality of the integral. To see that the integration constant is given by the local energy density, one starts from local energy conservation d 3 p (2π) 3 p μ u μ ( f − f iso ) = 0, and one rewrites the first term of this equation with the Landau matching condition u μ T ν μ = −ε u ν . The collision kernel for the first momentum moments F takes then the form − C[F] = −γ ε 1/4 −v μ u μ F + γ ε 5/4 −v μ u μ 3 .
To evaluate the kinetic equations (A21), (A22) for the harmonic coefficients F (n) , we need to expand the collision kernel C[F] to first order in the perturbation δ n . To this end, we note that the expansion (A20) defines a corresponding expansion T μν (t, r, θ) = T One can then determine the energy density ε = ε (0) + n =0 δ n e inθ ε (n) and the flow velocity by solving to linear order in δ n the Landau matching condition u μ T ν μ = −ε u ν . We proceed as follows: In the coordinate system (t, r, θ, z) with mainly positive metric g μν = diag −1, 1, r 2 , 1 and with the linearly independent unit vectors To determine the first order corrections ε (n) to the unperturbed local energy density ε (0) and the corrections u (n)r , u (n)θ to the radial and azimuthal component of the unperturbed flow vector (A34), we expand the Landau matching condition to first order in δ n , where T μν (n) (t, r, θ) = where v νR ν = − u−v ⊥ cos(φ r ) √ 1−u 2 and v νθ ν = v ⊥ sin φ r .
In this way, we have expressed the zeroth harmonics (A34), (A40), (A41), and the higher harmonics (A46), (A47) and (A48) in terms of integral moments of the solution F (n) and, a fortiori, the collision kernel is now given explicitly in terms of F.

Scaling of the initial conditions and of the equations of motion
For numerical evaluation, it is useful to write the initial conditions and equations of motion in dimensionless variables. To this end, one notes first that for very early times sufficiently close to τ 0 , F can be approximated by the free-streaming solution F(τ, x ⊥ ; φ, v z ) = 2δ(v z ) ε 0 τ 0 τ e − r 2 +τ 2 −2r τ cos φr For notational simplicity, this equation is written restricted to the second harmonic and to a gaussian initial profile, but explicit expressions can be given easily for all harmonics and for arbitrary radial profiles. One now defines the function F from F by rescaling all time and length scales with R, τ ≡ τ/R , x ⊥ ≡ x ⊥ /R and dividing out the prefactor R ε 0 τ 0 , To avoid overloading our notation, we drop the bars on dimensionless time and space coordinates immediately. The arguments of F will be understood in the following as being dimensionless, while the arguments of F are understood to be dimensionful. In particular, F(τ, x ⊥ ; φ, v z ) = 2δ(v z ) τ e −r 2 −τ 2 +2r τ cos φ r 1 + δ 2 [c (2) cos(2θ) + s (2) sin(2θ)] , where c (2) and s (2) are functions of the dimensionless r , τ and φ r , c (2) (r, φ r , τ ) = r 2 − 2r τ cos φ r + τ 2 cos(2φ r ), s (2) (r, φ r , τ ) = 2τ (r − τ cos φ r ) sin φ r .
When we insert (A52) into the equation of motion (3), we find Here, all space-time coordinates are understood as being dimensionless and the factors ε F on the right hand side are understood to be calculated from F. Since the energy density is defined as a momentum integral over F, we have ε = ε F ε 0 τ 0 R , where ε F is dimensionless. The prefactorγ combines the prefactor γ of the ITA in Eq. (A31) with a factor ε 0 τ 0 R 1/4 (arising from rescaling the energy density ε in the collision kernel (A31) and from changing from F to F), and with a factor R (arising from rescaling the spacetime derivatives and explicit 1 t -factors on the left hand side of (3)),