Scaling behavior of anisotropic flow harmonics in the far from equilibrium regime

We consider the development of anisotropic flow in an expanding system of particles undergoing very few rescatterings, using a kinetic-theoretical description with a nonlinear collision term. We derive the scaling behaviors of the harmonic coefficients vn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_n$$\end{document} with the initial-state eccentricities and the mean number of rescatterings, and argue that hexagonal flow v6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_6$$\end{document} should follow a nontrivial behavior, different from that of the lower harmonics. Our findings should be observable in experimental data for small systems.


Introduction
A highlight of the results from the ongoing experimental programs with heavy nuclei at the CERN LHC or the Brookhaven RHIC consists of the measurements of various azimuthal correlations between outgoing particles. In particular, the measured values of Fourier coefficients v n that quantify the anisotropies of the transverse emission pattern are interpreted as footprints of a strongly collective behavior [1], hinting at the creation of a medium close to local thermodynamic equilibrium.
According to the most widely accepted theoretical picture, the final-state anisotropies in momentum space reflect asymmetries of the spatial geometry of the "initial state"i.e. of the distribution of entropy density in the transverse plane after the nuclei have passed through each other [2,3] -, characterized for instance by "eccentricities" [4,5] n e in n ≡ − r n e inθ r n (1) a e-mail: borghini@physik.uni-bielefeld.de b e-mail: s.feld@physik.uni-bielefeld.de c e-mail: nkersting@physik.uni-bielefeld.de for n ≥ 2, and by generalizations thereof. 1 In this definition, r and θ are centered polar coordinates in the transverse plane, while the angular brackets denote an average over the entropy density. To be more specific, let us briefly quote a number of model findings, without any attempt at exhaustiveness (see also Ref. [6] for a compilation of experimental results from the LHC): • To a good approximation, the (integrated) elliptic flow v 2 and triangular flow v 3 scale linearly with the corresponding eccentricities [2][3][4][5]7,8]: 3 3 .
In collisions with a large impact parameter, i.e. a larger 2 , a cubic deviation to this behavior for v 2 has been reported [9].
• The quadrangular flow v 4 receives two kinds of contributions: a linear scaling with 4 and a quadratic dependence on 2 [5,7,[10][11][12]: Similarly, the pentagonal flow v 5 depends linearly on 5 and nonlinearly on the second-and third-harmonic eccentricities [5,11,12]: In both cases, the linear term is only visible in ultracentral collisions [8], while the nonlinear contribution dominates at larger impact parameters. • Starting with hexagonal flow v 6 , there appear more than one nonlinear terms at the "leading order" that mostly contributes in noncentral collisions [13][14][15]: v 6 K 6,6 6 + K 6,33 2 3 + K 6,24 2 4 + K 6,222 Note already that there are both quadratic and cubic contributions to v 6 , which will be of relevance for one of the findings of this paper.
The leading candidate theory for describing the evolution of the medium created in collisions of heavy nuclei at ultrarelativistic energies is nowadays relativistic fluid dynamics, and accordingly most of the studies quoted above were performed within that framework. The only exception is Ref. [8], which relies on a kinetic transport approach, however pushed into a regime where it "mimics" dissipative fluid dynamics. In the corresponding language, one finds that the linear response coefficients K n,n depend on the transport properties -mostly, the shear and bulk viscosity -of the expanding medium. This also holds for the nonlinear response coefficients K 4,22 , K 5,23 , K 6,222 , K 6,33 …, yet it was argued in Ref. [11] that they are less damped by viscous effects than the linear coefficients.
The collective-behavior picture underlying the fluiddynamical interpretation of the anisotropic flow data has however been challenged by the observations in the past few years of similar signals in collisions of "smaller systems" like proton-lead, deuteron-gold or even proton-proton, in events with a relatively large number of particles in the final state (see Ref. [16] for a recent review).
Accordingly, there has been a revival of models with "weak final-state collectivity" aiming at investigating generic behaviors of the anisotropic flow coefficients v n in a regime where the outgoing particles undergo in average very few rescatterings [17][18][19][20][21]. The present paper, which will remain at a semi-qualitative yet general level, is a further step in that direction, identifying a scaling behavior of the higher harmonics, in particular of the hexagonal flow v 6 , which to our knowledge has not been noted before. 2

Anisotropic flow far from equilibrium
Since our goal is to model a situation in which the final-state particles have only rescattered very little, a natural framework is to treat them as particles all along the system evolutionthey never form a continuous medium -and to resort to a kinetic theory framework, as we now detail.
Without loss of generality for our argumentation, we consider a single particle type and introduce its on-shell phase space distribution f (t, x, p), where three-vectors are denoted in boldface. The evolution of this distribution will generically be described by a kinetic equation of the form [23] with C[ f ] the collision term modeling the effect of rescatterings. 3 As initial condition for the evolution, we consider the following phase space distribution at some time t 0 , in which the dependences on the space and momentum variables factorize: where we again use polar coordinates (r, θ) for the projection of x onto the transverse plane, and discard the dependence on the longitudinal coordinate, which is irrelevant in the following. As hinted at by the notations, F(p) is a positionindependent momentum distribution, G(r ) depends only on the radial coordinate, while R is a length scale ensuring that 2 and˜ 3 are dimensionless. As noted by Teaney and Yan [4], the successive powers of r in the square brackets should actually be regulated by some cutoff function to ensure that the distribution f remains positive, yet we did not denote that regulator since it can remain unspecified for our purposes. Obviously, computing the eccentricities (1) with distribution (7) (instead of the corresponding entropy density, yet at the present stage this is only a matter of convention) gives n ∝˜ n for every n ≥ 2.
Letting now the initial distribution (7) evolve, we first consider the collisionless case, C[ f ] = 0. The corresponding solutions of the equation of motion (6) are free-streaming solutions with v the velocity corresponding to momentum p. The resulting anisotropic flow coefficients v n , 4 whose calculation involve integrations over the whole position space as well as over the momentum azimuthal angle φ p , are easily shown to be entirely determined by the initial transverse anisotropies of F(p), and in particular independent of the eccentricities n . Thus, if there is no anisotropic flow initially, as we shall from now on assume, there is none in the final state if the particles do not rescatter, as has long been known.
To turn on a small amount of rescatterings, we apply the idea whose various incarnations in the heavy-ion physics community went through the successive appellations "low density limit" [17,24], "far from equilibrium regime" [18] and more recently "eremitic expansion" [20] or "one-hit dynamics" [21], and write a solution of the kinetic equation (6) with collision term as where f (0) is the free-streaming solution (8) and f (1) a small correction term. Inserting this ansatz in Eq. (6) yields at once Multiplied by cos(nφ p ) and integrated over x and φ p , the term on the left hand side of this equation yields (up to a normalization factor) the negative of the time derivative ∂ t v n [18,22].
Integrating over time will then give the final v n . That is, specific integrals of the collision term C[ f (0) + f (1) ] yield the nth anisotropic flow harmonic v n .
Since we are interested in the development of the higher harmonics, and in particular in the nonlinear contributions, a natural choice for the collision term is Boltzmann's collision integral for elastic two-to-two scatterings, which we symbolically write in the form where the shorthand notation f ( j ) stands for f (t, x, p j ) while the terms w(i + j → k+l) involve transition probabilities and the necessary δ-distributions implementing energymomentum conservation. Note that in contrast to Ref. [18] we need not assume that the initial geometry is invariant under the x → −x transformation, nor that the involved interactions are parity non-violating. Given a model for the interaction, w will be proportional to some typical cross section σ . In turn, if one computes the total number of rescatterings taking place over the whole system evolution, it will also be approximately proportional to σ , as will be the average number of rescatterings per particleN resc. . The latter constitutes the dimensionless parameter that quantifies the smallness of f (1) relative to f (0) . When substituting f by f (0) + f (1) in this collision integral, we make use of the fact that f (1) is assumed to be a small correction and approximate [17,18,22] C f (0) + f (1) C f (0) .
That is, the free-streaming solution (8) fully determines the collision term and thereby the flow coefficients.
Performing the necessary calculations requires specific models for the as yet unspecified functions F(p) and G(r ) in the initial-state distribution (7) and for the interaction. General scalings can however already be predicted irrespective of any specific choice, which we now list.
• The multiplication of the isotropic term in one of the factors f (0) (i) in the integrand of the collision integral (11) with the term in¯ n cos[n(θ − n )] in the associated f (0) ( j ) yields a contribution to v n proportional to σ n , i.e. approximately proportional toN resc. n . With the values of the eccentricities relevant for heavy ion collisions, this is the dominant contribution to v 2 and v 3 , resulting in linear scalings of the form For n ≥ 4, other contributions to v n are likely to be as important, which we now discuss. • Besides the linear term in σ 4 , another contribution to v 4 is generated by multiplying together the terms in with a term quadratic in 2 yet linear in the mean number of rescatterings. Similarly, one finds v 5 ∼N resc. κ 5,5 5 +N resc. κ 5,23 2 3 , again with a contribution nonlinear in the initial-state eccentricities and linear inN resc. . • Coming to v 6 , one quickly sees that the approximation (12) will yield a contribution inN resc. To recover the latter, one must consider products f (0) (i) f (1) ( j ) in the integrand of the collision term, since f (1) contains a term in σ 2 2 cos(4θ) -which is reflected in the quadrangular flow (14). This will indeed yield a term in 3 2 , but the latter is also proportional to σ 2 : v 6 ∼N resc. κ 6,6 6 + κ 6,33 At very lowN resc. the linear and quadratic contributions will dominate, while the cubic term in the ellipticity 2 will only become meaningful for a larger number of rescatterings.
• Similarly, one can easily convince oneself that the setup consisting of the initial distribution (7), evolved with the Boltzmann equation (6) with collision term (11) does not generate any contribution to the harmonics v n≥7 involving only 2 and 3 at linear order inN resc. : the 2 2 3 contribution to v 7 is quadratic inN resc. ; v 8 will receive a term in 2 2 3 at orderN 2 resc. and a contribution 4 2 at orderN 3 resc. , and so on.
• Eventually, the same reasoning shows that the subleading contribution in 3 2 to elliptic flow v 2 , which is observed in fluid dynamical simulations, can also be recovered within our model, at orderN 2 resc. . More generally, the model predicts a contribution inN 2 resc. 3 n to v n for any n.
Summarizing our findings, we find that our model of a system of self-diffusing particles with an initially asymmetric transverse geometry is able to generate anisotropic flow coefficients v n with the same scaling dependence on the initialstate eccentricities as within a fluid-dynamical description, as seen from comparing Eqs. (13)- (16) with Eqs. (2)-(5). This behavior was already known in the literature [18,21,25].
What to our knowledge was never mentioned before regards the scaling of the generated flow harmonics with the rescattering cross section, or equivalently with the mean number of rescatterings per particle. 5 Thus, the "leading contributions", i.e. those stemming from the largest eccentricities 2 and 3 , to the successive Fourier coefficients v n scale with different powers ofN resc. . And, perhaps more interesting, starting with hexagonal flow v 6 there can be two or more such leading contributions to v n , which necessarily scale differently withN resc. : where the terms in 2 4 or 6 are assumed to be smaller. That is, the development of the contribution to v 6 from the ellipticity 2 necessitates more rescatterings than that of the triangularity 3 . Similarly, one easily finds v 8 4 2 , again assuming that the contributions to v 8 involving eccentricities p with p ≥ 4 are small.

Discussion
In this last section, we address two issues: first, are our results on the scalings with the average number of rescatterings N resc. , in particular that of Eq. (17), robust against changes of the setup which we considered? And second, is there a possibility to evidence these behaviors in experimental data?
If anisotropic flow is not present initially, but generated by rescatterings, then the corresponding harmonics will depend on the initial-state eccentricities n and onN resc. . Regarding the linear response of v n to n , we cannot think of a plausible scenario in which it would not already be generated at linear order inN resc. , i.e. v n = O(N resc. ) n , generalizing Eq. (13). Less straightforward are the nonlinear response behaviors, which we now discuss at length.
The emergence of scalings v n+ p ∝ n p at linear order inN resc. , as in Eqs. (14)- (15), is also straightforward in a description in which the collision term C[ f ] is at least quadratic in the single-particle distribution f , which is a natural feature in a picture in which the momentum anisotropies are generated by rescatterings of at least two partners.
The less trivial scaling behavior is that of Eq. (17), in particular the O(N 2 resc. ) dependence of the term in 3 2 . To investigate whether it is an artifact of our model or more general, we note that only two types of modifications are possible as long as one remains in a kinetic-theoretical framework with particle scatterings: changes in the initial-state distribution f (t 0 , x, p), which entirely determines the free-streaming solution f (0) (t, x, p), or of the collision term C[ f ]. In both cases, we want to see how a term in 3 2 might arise at first order inN resc. or somewhat equivalently in the interaction cross section σ , thereby invalidating Eq. (17).
Changing the initial-state distribution so as to spoil Eq. (17) is mathematically feasible, by assuming that F(p) contains a term in 2 . This would however mean that the initial momentum distribution already knows about the global geometry of the collision zone, which is problematic from the physics point of view. In turn, if the isotropic term G(r ) contains a term in 2 , then the latter will multiply both terms of Eq. (17), which is also not what we wish.
Accordingly, the only viable modifications to be considered are changes of the collision term C[ f ]. Sticking to the general structure of a collision integral involving f , 6 one quickly sees that the generalization of Boltzmann's ansatz (11) necessary to obtain a term 3 2 at first order inN resc. is to include a contribution of (at least) cubic order in f in the integrand. Two kinds of physical causes justify such contributions. On the one hand, one may include rescatterings with at least three particles in the initial state, in particular three-totwo scatterings, as can be found e.g. in Eq. (12) of Ref. [26]. Here, one should note that including two-to-three scatterings only would not help. In addition, three-to-two rescatterings will in fact generate the desired term in 3 2 at first order in the corresponding cross section σ 3→2 : whether the latter yields the leading contribution toN resc. , so that the term is indeed of order O(N resc. ) 3 2 , or elseN resc. is rather dominated by two-to-two rescatterings becomes a partly model-dependent issue.
On the other hand, even restricting oneself to elastic twoto-two rescatterings, one may still consider the generalization of the integrand accounting for Bose-Einstein enhancement or Pauli blocking, i.e. with factors of the form The inclusion of three-to-two scatterings and/or quantummechanical phase-space occupancy effects seems at first face to be relevant only if the initial state is that of a dense system. Naturally, when the latter expands, it becomes more dilute, and the terms beyond quadratic order in f in the collision integral become less important. Nevertheless, it is not clear to us whether a system created in high-energy collisions could be in a regime such that the emitted particles rescatter only very little, while at the same time being initially dense and interacting enough to lead to a breakdown of Eq. (17). "Small systems" are the most likely candidates for such a departure, provided the initial density is big. In any case, the relative importance of more-than-quadratic terms and their influence on the scaling behavior could be tested in numerical simulations with transport codes in which they can be switched on or off at will, like 2 ↔ 3 scatterings in BAMPS [26] or quantum effects in other codes [27].
Let us now discuss where in experimental data the scaling behaviors (13)- (17) could possibly be at play and measurable.
Surprising though it may seem, let us first deal with larger systems, in which fluid dynamics is routinely applied to describe the evolution. On the one hand, the single-collision regime might be applicable to particles in given regions of phase space, e.g. at high transverse momentum [20], or to specific particle types, like bottomonia -for which the picture is rather a negative one: a collision means destruction. On the other hand, Eqs. (13)-(17) may also be relevant in phenomenological analyses of the bulk of particles. More accurately, these behaviors play a role in the pre-hydrodynamized stage, which in modern hybrid descriptions is often modeled by a transport cascade [28]. Indeed, this short kinetic period, which only involves rather few rescatterings, will transform "pre-early transport eccentricities", taken from a model for initial conditions, into some early anisotropic flow, which becomes part of the initial condition for the fluid-dynamical stage of the hybrid description. Following Eq. (17), the early generated v 6 will suffer from a deficit in second-order eccentricity 2 , which will be propagated by the subsequent evolution until the final state. This could then affect attempts at evidencing the scaling (5) and at interpreting it within a purely fluid-dynamical framework, since the coefficients K 6,33 and K 6,222 will contain a pre-hydrodynamization component, whose relative size possibly varies across centralities. This possibility is a further incentive to investigate the scaling behaviors (13)- (17) in transport models.
Eventually, the natural place where Eqs. (13)- (17) are to be looked for is in small systems, in which the applicability of fluid dynamics is most questionable. The biggest issue is of course that the anisotropic flow coefficients in such systems are small. We believe that the more trivial scaling behaviors (13)-(15) should be rather "easily" observable. Note in particular that the unknown mean number of rescatterings cancels in the ratio of two different harmonics v 2 -v 5 , so that one can separate the influences of eccentricities andN resc. , where one can expect that the latter should scale like the cubic root of the charged particle multiplicity (dN ch. /dη) 1/3 . In turn, we are aware that in small systems v 6 will be at the border of what is measurable with reasonable uncertainties, so that whether measurements allowing to test Eq. (17) are feasible is not warranted. Nevertheless, we think that confirming the scaling behavior (17) would yield further confidence in the determined value ofN resc. , at the same time evidencing a nice instance of nonlinearity. Conversely, as we have already discussed above, departure from that behavior might hint at a dense initial state, possibly saturated, which is certainly not an uninteresting result and is worth investigating.