On differences between even and odd anisotropic-flow harmonics in non-equilibrated systems

To assess how anisotropic transverse flow is created in a system out of equilibrium, we compare several kinetic-theoretical models in the few-rescatterings regime. We compare the flow harmonics 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} from three types of transport simulations, with either 2→2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\rightarrow 2$$\end{document} or 2→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\rightarrow 0$$\end{document} collision kernels and in the former case allowing the particles to rescatter several times or not, and from analytical calculations neglecting the gain term of the Boltzmann equation. We find that the even flow harmonics are similar in all approaches, while the odd ones differ significantly. This suggests that while even 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} harmonics may to a large extent be due to the anisotropic escape probability of particles, this is not the predominant mechanism underlying the odd 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} coefficients.


I. INTRODUCTION
The charged hadrons produced in collisions of heavy nuclei at high energy show a characteristic azimuthally asymmetric transverse emission pattern [1].This anisotropic flow, usually quantified in terms of coefficients in the Fourier expansion of the transverse momentum distributions [2], has also been observed in so-called smaller systems, namely proton-and deuteron-nucleus or even proton-proton collisions with large multiplicities [3].
The anisotropic flow results have been essential for singling out relativistic hydrodynamics [4,5] as the model of choice for describing the dynamics of the system created in heavy-ion collisions, which is then modeled as a continuous medium, whose initial geometrical asymmetry is converted in the evolution into the final state momentum anisotropy [6].Indeed, relativistic fluid dynamics can describe satisfactorily a large amount of data on anisotropic flow and related azimuthal correlations [1,[7][8][9].
However, the question of the applicability of fluid dynamics is still under discussion, in particular when the number of emitted particles is small [10,11].Thus, alternative descriptions based on microscopic kinetic transport theory, which is known to reproduce fluid-dynamical results when particles undergo many rescatterings [12], are being explored again, in particular with a view to small systems.A number of recent attempts start from semi-realistic initial geometries, which allow to isolate the flow harmonics of interest and study their origin [13][14][15][16][17][18][19][20].
In one of the more realistic transport studies [21], within the AMPT approach, it was claimed that "the majority" of the measured anisotropic flow signal (for elliptic flow v 2 and triangular flow v 3 in Au-Au collisions at RHIC energy) is not due to the numerous rescatterings of the produced particles, but could rather be dominated by those particles that escape the asymmetric system geometry without having scattered.Primitive versions of this "escape mechanism" scenario had been considered earlier with simple initial states allowing analytical calculations with simplifying assumptions [22,23], and also used for an early estimate of the v 2 of J/ψ quarkonia [24].
Yet recent findings cast some doubt on the efficiency of the anisotropic-escape picture in the regime of very few rescatterings, especially regarding v 3 .Thus, it was found in Ref. [18] that the value of v 3 -to be accurate, of energy-weighted triangular flow -in kinetic models at low opacity, i.e. when particles rescatter very little, depends significantly on the collision kernel of the Boltzmann equation: triangular flow (divided by the initial triangularity) comes out negative in an effective kinetic theory of QCD, while it is positive in the relaxation-time approximation.In contrast, the behavior of elliptic flow seems to be more robust across scenarios.
In this paper, we want to further explore the production of anisotropic flow in the regime of very few rescatterings, in particular with a view to testing the anisotropic-escape scenario.For that purpose, we employ numerical transport simulations with various collision kernels, in particular with elastic binary collisions (Sect.II), complemented with analytical calculations that only account for the loss term of the Boltzmann equation.We then compare in Sect.III the results for the v n coefficients in our various approaches and with those of the recent literature, before concluding in Sect.IV.
Since we focus on systems with very few rescatterings, the flow coefficients are at times very small.Accordingly, their values in transport simulations are likely to be affected by numerical fluctuations.High statistics are needed to counteract this noise, which is why we restrict ourselves to a two-dimensional system, to keep the computing time in reasonable bounds.This restriction will be further examined in Sect.IV.

II. METHODS
In order to investigate the importance of the "escape mechanism" for the production of anisotropic flow when particles undergo very few rescatterings, we perform four types of calculations.On the one hand, numerical simulations with a transport code, with two different colli-sion kernels: first a 2 → 2 kernel that implements elastic binary collisions, and gives as reference the "total" anisotropic flow produced in a semi-realistic system.Secondly, a "single-hit" version using the 2 → 2 kernel but in which particles that have already scattered once are no longer allowed to interact.Thirdly, a 2 → 0 collision kernel, such that the resulting flow is that of the particles that escaped the system without scattering.On the other hand, we perform analytical calculations within kinetic theory, using only the loss term of the collision kernel in the Boltzmann equation, and working at linear order in the cross section: this provides a controlled approximation to the 2 → 0 scenario, which itself includes all orders in the cross section.
We begin with introducing the analytical approach (Sect.II A), together with the initial conditions we use for both analytical and numerical calculations.We then briefly present our transport setups (Sect.II B).All calculations are performed with massless identical (yet distinguishable) particles, which propagate in two dimensions only, corresponding to the transverse plane in a high-energy nuclear collision.Two-dimensional vectors are denoted in boldface.Throughout the paper we use the convention h = c = 1, and (r, θ) denote polar coordinates in the transverse plane, with their origin at the center of the system in its initial state.

A. Analytical approach
In our analytical calculations, we characterize the particle system by a classical on-shell phase space distribution f , which obeys the relativistic Boltzmann equation Instead of a full collision kernel with detailed balance, we shall only consider the loss term of binary scatterings with E p the energy of the particle with momentum p, v rel. the Møller velocity, and σ the total cross section.Note that this implies that energy, momentum and particle number are not conserved in the evolution.For massless particles in two dimensions, v rel.= 1 − cos(ϕ p − ϕ 1 ) where ϕ p resp.ϕ 1 is the azimuthal angle of momentum p resp.p 1 .
The "observables" we study are the Fourier coefficients quantifying anisotropic flow [2], in particular their time evolution.In terms of the phase space distribution f , the momentum-integrated coefficients are given by where the denominator is simply the total number of particles N (t) at time t.Differentiating this expression with respect to time gives two contributions, from the derivatives of the numerator and denominator respectively: Using the Boltzmann equation to replace ∂ t f in the integrand, the term involving the spatial gradient of f gives zero after integrating over x, since f vanishes at infinity.There remains only the contribution from the collision term, which at leading order is a priori linear in σ.If we restrict ourselves to this linear order, as we do from now on, then we may neglect the change in N (t) induced by the (particle-number non-conserving) rescatterings in the denominator in the first line of Eq. ( 4), i.e. approximate N (t) N (0), which we shall more briefly denote by N .In addition, we may also neglect the evolution of the phase-space density induced by rescatterings in the integrand of the collision term in the numerator.That is, we replace f (t, x, p) by the free-streaming distribution f f.s.(t, x, p) that coincides with f in the initial state [13,22,23]: where v ≡ p/|p| while f (0) (x, p) denotes the initial distribution (at t = 0), to which we come back hereafter.In the second line of Eq. ( 4), ∂ t N (t) is of order O(σ) (or higher).In absence of initial anisotropic flow in the system, v n (t) is also of order O(σ), so that the whole term is at least quadratic in σ: accordingly, we shall neglect it hereafter.Note however that this term contributes at linear order in σ, and thus may not be dropped, if there is some anisotropic flow in the initial state.
All in all, we replace the evolution equation ( 4) with valid to linear order in σ, irrespective of the choice of collision term -as long as the latter is O(σ).Inserting the loss term (2) as collision kernel and integrating over time yields where the last line defines the angle-averaged local production rate of v n [18], which we shall discuss in Sect.III E. Note that in these expressions we explicitly assumed v n (t = 0) = 0 in the initial state.
In our analytical approach the flow coefficients (7) depend directly on the initial phase space distribution f (0) via Eq.( 5).Let us now discuss our choice for the latter, both for the analytical calculations and the numerical simulations.First, we assume that the initial phase space distribution factorizes into the product of the particle number density, which determines the geometry, and a position-independent momentum distribution: where we assume that G is normalized to unity when integrating over the whole two-dimensional momentum space.This factorization assumption makes our analytical calculations tractable, and enables us to derive analytical formulas for the flow coefficients for the geometrical profile (9).As we shall discuss again in the following, the assumption is however not innocuous, especially for the odd flow harmonics.We take G to be isotropic in momentum space, to ensure the absence of initial anisotropic flow.Departure from this assumption can be accounted for rather easily, by introducing a Fourier expansion of G(p) [25], but leads to lengthier expressions for the flow coefficients -whose evolution at linear order in σ is no longer governed by Eq. ( 6) as mentioned above.
In position space, we choose as initial density a distorted Gaussian distribution1 with N the number of particles and R the typical system size, in units of which we shall measure lengths or time.This form or closely related ones was used extensively in recent studies [14][15][16][17][18][19]26], as it allows one to introduce at will in the initial state different and independent types of "eccentricities" [27][28][29] ε n e inΦn ≡ − r n e inθ r n , where the angular brackets stand for an average over the transverse plane with some weight, which in the present paper will be the particle-number density.Equation ( 9) yields at once Φ n = 0 -which we may assume without loss of generality since we shall always consider only a single non-zero ε n at a time -and that is for the first harmonics ε 2 = ε2 /4, ε 3 = ε3 / √ 2π, ε 4 = 3 ε4 /4, and so on.Note that the parameters εn should not be too large, to ensure that the phase space distribution remains non-negative: typically, in case only a single eccentricity is considered, εn should be such that ε n remains smaller than ε n,max 0.35.In our calculations, both analytical and numerical, we choose εn such that ε n = 0.15 or smaller.

B. Numerical simulations
For our simulations with elastic binary rescatterings, we use the same implementation of the two-dimensional covariant transport algorithm of Ref. [12] as in Ref. [30], to which we refer for further details.Here we just recall that the N massless particles are modeled as N p Lorentzcontracted hard spheres -or rather hard disks, since they are two-dimensional -with radius (N/N p )σ/2, where σ is the total cross section of the "physical" particles.Collisions between test particles are determined by a geometric criterion and the scattering angle is deterministic.N p and σ are always chosen such that the system remains dilute enough, i.e. the mean inter-particle distance is at least one order of magnitude smaller than the mean free path mfp .
For the simulations with the 2 → 0 collision kernel we use the same transport algorithm as in the 2 → 2 case with small modifications.We introduce labels "active" and "inactive" for each test particle, such that a collision can only take place between two "active" particles, after which they become "inactive" and are no longer propagated for the remainder of the simulation.Eventually, observables like the anisotropic flow coefficients are determined with the "active" particles only.
An important difference between this 2 → 0 model and the analytical approach is that the phase-space distribution in the simulations is affected by rescatterings, i.e. the transport simulations include all orders in the cross section.Thus, we may depart from the few-collision regime in the simulations and investigate what happens when most of the particles disappear due to rescatterings.
Eventually, we also consider a third variant, which we shall refer to as "single hit" model, in which particles scatter with the 2 → 2 kernel, but may undergo at most one collision.That is, after their first rescattering -and the corresponding change in the momenta of the two participants -, particles become "transparent" and stream freely through the system.The difference with the 2 → 0 scenario is that all particles are now taken into account when computing anisotropic flow, irrespective of whether they have undergone zero or one collision.
In Ref. [21] the authors used a similar approach with 2 → 2 and 2 → 0 collisions.The difference to our 2 → 0 model is that in their study, particles that underwent a collision are still "active", but after each collision their momentum azimuths are randomized.Thus, these particles do indirectly contribute to the generation of anisotropic flow in the azimuth-randomized version of AMPT [21]. 2 → 2 2 → 0 loss term FIG. 1. Mean number of rescatterings per particle over the system evolution as a function of the inverse Knudsen number estimated in the initial state, Eq. ( 12), for the 2 → 2 (green squares, fit with Nresc.≈ 0.529 Kn −1 ) and 2 → 0 (blue circles, fitted with a quadratic ansatz: dashed line) scenarios.The red line Nresc.= Kn −1 / √ π is the prediction of the analytical approach.
A crucial ingredient for the comparison with our analytical calculations is the preparation of the initial state of the numerical simulations.The test particle positions are sampled from the distribution function ( 9), while for their momenta we use a Boltzmann distribution with a position-independent temperature -in contrast to Ref. [20].Since the simulations are performed with a finite test particle number N p ranging between 2 × 10 5 and 2 × 10 6 , neither perfect isotropy in momentum space nor uniformity of the momentum distribution across the whole geometry can be achieved.To improve the situation, for each initial geometry we perform N iter.iterations in which the particles keep the same position but with a different realization of the momentum distribution.The results we present are averaged over these iterations, which is expected to diminish fluctuations by a factor √ N iter. .Since the simulation time grows with N 3/2 p , performing multiple iterations with less test particles is computationally less costly than performing a single simulation with N iter.N p particles. 2  Starting from Eq. ( 9), the average particle-number density per unit surface is N/4πR 2 .Using the latter to define a mean free path mfp , we quantify the rarity or abundance of rescatterings by the Knudsen number with the help of which we shall express the equations resulting from the analytical calculations.In contrast, the results of numerical simulations will be presented not at 2 In our simulations, N iter.Np is always larger than 10 9 .2 → 0 2 → 2 loss term FIG. 2. Time dependence of the cumulative number of rescatterings per particle for systems with in total Nresc.≈ 0.14 (full) or 0.08 (dashed) at large times, for simulations with the 2 → 2 (green) and 2 → 0 (blue) models, and using Eq. ( 13) (red).
fixed Kn, but rather at fixed mean number of rescatterings per particle N resc.over the whole evolution -in practice, until t/R = 30.We shall mostly present results for N resc.≈ 0.02, well in the few-rescatterings regime, and 0.14 -for which the approximation becomes less justified -, as well as N resc.≈ 0.35 in Appendix C.
In the 2 → 2 scenario, N resc.nicely scales with Kn −1 , see Fig. 1.Note that there are slightly less (about 8%) rescatterings in our simulations than what would be expected analytically.This is due to the finite time step of the transport code, and to the fact that a given particle is allowed to scatter only once per time step, so that we miss collisions,3 mostly in the densest regions of the system.That is, the effective Kn −1 in the simulations is actually smaller than that computed from the input parameters, which is a first motivation for presenting numerical results in terms of N resc.instead.A second reason for using the mean number of rescatterings per particle is that it turns out that it is the correct scaling variable for comparing systems in the 2 → 0 scenario, as will be discussed hereafter in Sect.III A and III B.
Nevertheless, it is clear that a given N resc.requires a larger Kn −1 , i.e. cross section, in the 2 → 0 and single-hit models than in the 2 → 2 simulations, since particles can never scatter twice in those scenarios.This in turn means that the collisions tend to occur earlier in the 2 → 0 and single-hit simulations than in the 2 → 2 model, as is illustrated in Fig. 2 for azimuthally symmetric systems with in total N resc.≈ 0.08 (dashed) or 0.14 (full lines) collisions per particle.Accordingly, the geometry of the system at the time of the rescatterings varies across the setups.For instance, since the initial asymmetries in the geometry relax as the system expands, one may expect that at the time when anisotropic flow develops -say roughly for t/R ≤ 2 -the system is somewhat more isotropic in the 2 → 2 simulations than in the other ones, which impacts the anisotropic flow coefficients.Anticipating on our findings, this effect does not seem to play a major role.
In Fig. 2 we also show the time dependence of the number of rescatterings within the analytical approach of Sect.II A, i.e. using the free-streaming phase-space distribution all along the evolution.For the initial distribution ( 9) with vanishing eccentricities one finds ) with I 0 and I 1 modified Bessel functions of the first kind.Choosing the value of Kn −1 such that it yields the same final N resc.as in the numerical simulations, we see that this formula gives an extremely good approximation to the results in the 2 → 2 model.

III. RESULTS
In this Section we present our results for the flow harmonics v 2 , v 3 , v 4 , and v 6 for systems with the initial geometry (9).Calculations with a slightly different initial profile, whose results are in qualitative agreement with the findings of this Section, are given in Appendix D.

A. Elliptic flow
Let us start with elliptic flow v 2 [6].As initial geometry we consider the profile (9) with all εj = 0 except for ε2 , chosen such that ε 2 = 0.15 (up to numerical fluctuations in the simulations).The time dependence of v 2 in the transport approach is shown in Fig. 3 within the 2 → 2 (green), 2 → 0 (blue) and single-hit (orange) models, for N resc ≈ 0.02 (left panel) and N resc ≈ 0.14 (right panel). 4 At t = 0 we indicate as an error bar the typical value 1/ 2N iter.N p of v 2 induced by numerical fluctuations in the initial state.We also show in red the result from the analytical calculation, namely where the value of Kn −1 is chosen such that it gives the same N resc.as in the numerical calculations.Note that 4 Results in systems with Nresc.≈ 0.35 are shown in Fig. 12.
To quantify the deviation between the various approaches, we fitted our results from transport simulations, shifted to v 2 (t = 0) = 0 for a better comparison, with respective Padé approximants to wash out the numerical fluctuations, especially at early times.A drawback from the approximation is that the fits are dominated by the values for t/R > ∼ 1, so that the early time behaviors are not necessarily captured correctly.Using these fits, we computed the ratios of the v 2 values in the 2 → 0 scenario either to those of the 2 → 2 and single-hit models or to the analytical value (14), and show these ratios in the narrow lower panels in Fig. 3.
The profiles of v 2 (t) are similar in the four approaches, with a slow onset, followed by an almost linear rise, that eventually saturates.v 2 reaches its maximum value for t/R ≈ 2, and decreases a little afterwards, barely in the 2 → 2 and single-hit scenarios.Remarkably, the overall shape of v 2 (t) is the same for the small numbers of rescatterings considered here as in the fluid-dynamical limit, illustrated e.g. in Ref. [31] (Fig. 3, with a slightly different geometry).
More importantly for the purpose of this paper, the elliptic flow built up in the 2 → 0 model differs at most by 20% from that in the "full" 2 → 2 case.In addition, the results of the 2 → 0 scenario are extremely well reproduced by the analytical calculations accounting only for the loss term for N resc.≈ 0.02.The agreement is less impressive but still very good at the larger N resc., which is easily understood: The analytical results are derived at linear order in the cross section, or equivalently Kn −1 .As the latter increases, higher order contributions to v 2 , which are always present in the 2 → 0 simulations, become more sizable, and lead to the departure between the analytical results and the 2 → 0 values.Indeed, we have shown in Ref. [20] -yet only for early timesthat pushing the analytical calculation to higher order in σ improves the agreement with the 2 → 0 results.In contrast, for t/R ≥ 2, when fewer collisions take place, the results of both approaches are again very parallel.
Eventually, the results from the single-hit scenario show a non-systematic trend with varying N resc. .When the number of rescatterings is very small (N resc.≈ 0.02), the single-hit v 2 (t) is intermediate between the 2 → 0 and 2 → 2 results.This seems consistent with the intuition that the single-hit model captures part of the gain term of the Boltzmann equation -since colliding particles are redistributed in momentum space -, but not the whole of it, as particles can scatter at most once.However, when N resc.increases, the single-hit results for v 2 (t) depart more strongly from those of the 2 → 2 cascade, and they are now further away from them as those from the 2 → 0 scenario, see right panel of Fig. 3. FIG. 3. Top: Time dependence of elliptic flow v2 in systems with initially ε2 = 0.15 and on average Nresc.≈ 0.02 (left) or 0.14 (right) rescatterings per particle.The green curves are for systems with elastic binary scatterings, the blue lines for the 2 → 0 scenario, the orange curves for the single-hit model, and the red ones show the analytical result (14).The bottom panels show ratios of the curves from the upper panels.FIG. 4. Ratio v2(t)/Nresc.in systems with initially ε2 = 0.15 and on average Nresc.≈ 0.02 (full lines) or 0.14 (dashed) rescatterings per particle for the three scenarios of the transport cascade: 2 → 2 (green), 2 → 0 (blue), single hit (orange).
The somewhat different behavior of the single-hit model, for which we could not find an easy explanation, is also illustrated in Fig. 4, which displays v 2 (t) scaled by the total number of rescatterings for the three scenarios of our transport code and for the two values N resc.≈ 0.02 and 0.14.This figure shows that to a very good approximation v 2 (t) ∝ N resc.holds in the full 2 → 2 simulations --it is then equivalent to v 2 ∝ Kn −1 , see Fig. 1, i.e. v 2 ∝ σ -and the 2 → 0 model, but the scaling is less good, although still satisfactory, for the single-hit case.In Fig. 11 in Appendix B we show for the sake of completeness the ratio of v 2 (t)/Kn −1 , i.e. essentially of elliptic flow over the cross section, for the same simulations as in Fig. 4. For the 2 → 0 and single-hit scenarios, the curves corresponding to systems with N resc.≈ 0.02 and 0.14 are far apart from each other, which shows that N resc. is indeed a better scaling variable than the inverse Knudsen number Kn −1 for those simulations (at least as far as anisotropic flow is concerned).
All in all, we find that in the few-rescatterings regime most of the v 2 signal may be ascribed to the processes modeled by the loss term of the Boltzmann equation.That is, the elliptic flow in the final state seems to arise to a large extent from the anisotropic survival probability of the particles as they propagate through the system [22,23], as advocated in the "escape mechanism" picture [21].

B. Triangular flow
We turn next to triangular flow v 3 [27], using now an initial geometrical profile (9) with only a non-zero ε3 , such that ε 3 = 0.15.The results of our various calculations for the time dependence of v 3 are displayed in Fig. 5, for systems with N resc.≈ 0.02 (top left), 0.08 (top right) or 0.14 (bottom).
A first striking feature is that v 3 identically vanishes in the analytical approach if it is zero initially.As we show in Appendix A and discuss again in Sect.III E, this is due to a cancellation between different regions in the special case -which we consider throughout the paper -where the local momentum distribution is the same at every point of the transverse plane in the initial condition.To be more precise, one finds that v 3 , and more generally every odd flow harmonic, is zero at first order in σ, but at higher orders it can be non-zero [20].
As to the results of transport simulations, we see a number of differences with those for elliptic flow.First, the v 3 signal is an order of magnitude smaller than v 2 , so that the curves are more affected by the numerical fluctuations, in particular in the initial state. 5Secondly, the simulations within the 2 → 0 model give a clear nonzero signal, in contrast to the analytical result.This hints that in the 2 → 0 simulations, which include all orders in the cross section, v 3 arises at a higher order in σ.
Thirdly, the results of the 2 → 0 scenario clearly do not resemble those of the 2 → 2 model.For N resc.≈ 0.02, the 2 → 0 results lie about a factor 1.5 below, while they are larger for N resc.≈ 0.08 and 0.14.
Eventually, the results of the single-hit model for v 3 (t) again show no clear trend in comparison to the other two numerical models.At N resc.≈ 0.02 they closely resemble the results of the 2 → 2 computations -the overshooting is probably due to the initial noise.But at larger number of rescatterings they are closer to the outcome of the 2 → 0 simulations, which makes it difficult to draw any conclusion.
In Fig. 6 we compare systems with different number of rescatterings by scaling v 3 by N resc.(left) or N 2 resc.
we discard them from the comparison and only look at N resc.≈ 0.08 and 0.14.For the 2 → 0 scenario, the plots hint at a scaling behavior v 3 ∝ N 2 resc., different from that found for elliptic flow.Regarding the 2 → 2 (and even more the single-hit models), the plots are rather inconclusive, and both scalings with N resc.and N 2 resc.seem almost acceptable.Let us note that studies focusing on the final value of v 3 , at the end of the evolution, have found v 3 ∝ N resc.(or equivalently v 3 ∝ Kn −1 ) at small N resc. in systems with elastic binary scalings [20,31].
Several recent studies investigated the "energy weighted triangular flow" v E 3 , i.e. the third Fourier coefficient of the transverse energy distribution, instead of the particle-number weighted coefficients [17][18][19].As shown in Fig. 7, v E 3 -computed in the same systems as used for Fig. 5 -again differs a lot in the 2 → 2 and 2 → 0 scenarios.This is especially true at times t > ∼ R. In turn, the single-hit results are quite close to the 2 → 2 values at N resc.≈ 0.02, but at higher N resc.they tend to be systematically larger.At earlier times t < ∼ R, the results with the three scenarios are more similar, but this is possibly a coincidence, as part of that early behavior is driven by numerical fluctuations: due to the finite number of particles, it is impossible to impose that the momentum distribution be exactly isotropic and identical everywhere in the transverse plane, so that the numerical realizations differ from the idealized setup.
In summary, and in strong contrast to the findings of Sect.III A, we find that for v 3 the 2 → 0 scenario differs significantly from the 2 → 2 model.In parallel, the triangular flow from the analytical approach considering only the loss term at first order in σ is also at variance with the results of numerical simulations. 6This is a strong hint that the final state of the individual rescatterings, modeled by the gain term of the Boltzmann equation, plays a crucial role: That is, the observed v 3 is not carried predominantly by particles that underwent no rescattering and escaped anisotropically from the medium.

C. Quadrangular flow
With quadrangular flow v 4 , the situation is again simpler than for v 3 .Anticipating on what we shall now present, the overall trend is the same as for elliptic flow v 2 : the results of the numerical 2 → 2, 2 → 0 and single-hit simulations and those of the analytical approach nicely agree when the number of rescatterings is (very) small, hinting at the dominant role of the escape mechanism for v 4 in this regime.
Starting with v 4 , a new possibility appears, namely that the produced anisotropic flow harmonic v n can arise not only because of the spatial harmonic ε n , but also due to nonlinear effects mixing other eccentricities. 7Thus, v 4 may be caused not only by the "quadrangularity" ε 4 , but also by the "ellipticity" ε 2 [14,29,[33][34][35].
Indeed, our analytical calculation for v 4 assuming only with initially a thermal momentum distribution with a positiondependent temperature -, the results of analytical calculations for v 3 (t) at order σ 2 but restricted to early times t ≤ R are found to be of the same magnitude as those of numerical computations, but the shape (which is affected by numerical noise) is not reproduced. 7To be more accurate, according to our present knowledge the lower flow harmonics vn with n ≤ 3 are only minimally affected by such nonlinear effects involving eccentricities ε k with k = n.a non-vanishing ε4 in the initial state yields Assuming instead that only a non-vanishing ε2 is initially present, we obtain Obviously, the terms on the right-hand sides of these equations add up if the initial state contains both ε2 and ε4 .These analytical results are compared to those of numerical simulations with both 2 → 2 (green) and 2 → 0 (blue) collision kernels in Fig. 8: the plots in the top panels are with ε4 = 0, such that ε 4 = 0.15, and all other εk = 0, while the bottom panels -in which we also show the results from simulations in the single-hit scenario (orange) -are for a non-zero ε2 (with ε 2 = 0.15) and vanishing other eccentricities. 8Figure 8 displays the time evolution of v 4 for systems with N resc.≈ 0.02 (left) or 0.14 (right) rescatterings per particle, while results for N resc.≈ 0.35 are shown in Fig. 13.Overall, the results in Fig. 8 show that in the case of 8 These simulations with only an ellipticity ε 2 = 0.15 are actually the same as used for v 2 in Sect.III A. quadrangular flow v 4 , either from ε 2 or from ε 4 , the 2 → 0 model represents a very good approximation of the 2 → 2 collision kernel for low N resc. .In turn, the nice agreement with the analytical results reinforces that statement and shows that v 4 is proportional to σ in that regime.Indeed, the less good agreement of the "loss term" results with the 2 → 0 kernel for N resc.= 0.14 can be attributed to the limitation of the analytical calculations to linear order in the cross section.At both values of N resc.and for collisions with an initial ellipticity, the v 4 values from the single-hit model also roughly match those of the 2 → 2 simulations, although less so at the larger N resc. .Although the results of Fig. 8 suggest that v 4 behaves as v 2 , in that it seems to be mostly driven by the particles that did not collide -at least in the low N resc.regime -, still there are important differences.A first one, to which we shall come back in Sect.III E, is that v 4 changes sign over time, while v 2 does not.A second difference is that while the overall shape of v 2 (t) is roughly the same in the few-rescatterings regime and in the fluid-dynamical limit, this does not hold true for v 4 (t).Indeed, we find that for a larger number of rescatterings (N resc.> ∼ 5, with the 2 → 2 collision kernel, since the 2 → 0 scenario makes no sense in that case) the v 4 resulting from an initial ε 4 > 0 is positive at late times, as found also e.g. in Refs.[17,31], 9 but contrary to the behavior of the upper panels of Fig. 8.This means that the linear scaling with N resc. of the "final" v 4 observed in Fig. 8 breaks down at larger cross sections.Note that a negative v 4 -more accurately, v E 4 -for a positive ε 4 in the few-rescatterings 9. Time dependence of hexagonal flow v6 in systems with on average Nresc.≈ 0.02 per particle and with different initial geometrical profiles, as described in the text.
regime was also found in Ref. [19], yet with a different collision kernel based on the relaxation time approximation.This difference in the collision kernel may explain why we do not find the same behavior at early timesnamely a negative v 4 -in case the system is initially deformed elliptically (ε 2 = 0, ε 4 = 0).
All in all, it seems that in the few-rescatterings regime v 4 , either resulting "linearly" from an initial quadrangularity ε 4 or nonlinearly from an initial ellipticity ε 2 , behaves like elliptic flow v 2 , i.e. it largely arises from the anisotropic escape of particles.Interestingly, the contributions from ε 2 and ε 4 to v 4 are of the same order of magnitude, and in the small N resc.regime they are of opposite signs.Accordingly, the two contributions can partly cancel each other and lead to a v 4 value at large times that can lie in a wide range of values.In particular, it is possible to obtain a negative v 4 value.

D. Hexagonal flow
Going beyond v 4 , we can guess qualitatively in analogy to our study of v 3 what we would find for v 5 : since it is an odd harmonic, the analytical approach gives zero at linear order in σ.In turn, this means that in the 2 → 0 scenario v 5 arises at order N 2 resc., while it is proportional to N resc. in the 2 → 2 model, so that we would find discrepancies between the two types of transport simulations.We did not attempt to perform such simulations, which would require new sets of events with the appropriate controlled initial geometry.
Here we present in Fig. 9 results for v 6 , which is at the limit of what we can do numerically with reasonable control on the signal when N resc.≈ 0.02, while exploiting simulations that were already used for v 2 or v 3 .Indeed, an interesting feature of v 6 is that it can result from different initial geometries [36][37][38], in particular with only a hexagonal deformation (linear response v 6 ∝ ε 6 , dotdashed line), only an initial triangularity (quadratic response v 6 ∝ ε 2 3 , dashed lines), only an initial ε 2 (cubic response v 6 ∝ ε 3  2 , dotted lines), or with both initial ε 2 and ε 4 (mixed quadratic response v 6 ∝ ε 2 ε 4 , full line).In every setup the only non-zero ε n are set to 0.15.The numerical results with an initial ε 2 resp.ε 3 are from the same simulations as in Sect.III A resp.III B. We did not attempt to perform simulations with an initial non-zero ε 6 nor with both ε 2 and ε 4 (and aligned symmetry planes Φ 2 and Φ 4 , as assumed for the analytical curve).
Similarly to what we found for v 2 and v 4 , the results for v 6 stemming from an initial ε 3 agree rather well across the three scenarios of this paper in the fewrescatterings regime.This agreement should be contrasted with Sect.III B, in which the same initial setup yielded very disparate results for v 3 .This reinforces our main conclusion of the paper regarding the different "origins" of the even and odd flow harmonics.
As regards the v 6 from an initial ε 2 , the results from numerical simulations are extremely small but seem to be non-zero and consistent in the 2 → 2 and 2 → 0 models.In contrast, the analytical results in that case are exactly zero: as was pointed out in Ref. [14], in a model with only binary collisions and no quantum-statistical effects, a contribution in ε 3  2 to v 6 (or to v 2 ) can only arise at order σ 2 , not at linear order in σ as considered here.
Eventually, the analytical results for initial geometries with either ε 6 = 0.15 or ε 2 = ε 4 = 0.15 are of the same typical magnitude as those for ε 3 = 0.15.As in the case of v 4 , the signal changes sign (here twice) over the system evolution.

E. Local production rate of anisotropic flow
To probe the temporal and spatial origin of the anisotropic flow buildup better, we study the production rate of each flow harmonic as a function of time and position [18,19].This local production rate of v n , averaged over the polar angle of the production point, is quantified by D n (t, r) introduced in Eq. ( 7), from which its expression can be read off. Figure 10 shows the results of our analytical approach for D 2 , D 3 , and D 4 , for the setups of Sects.III A-III C with N resc.= 0.02.Note that we display D n (t, r) multiplied by r, so that the production rate ∂ t v n (t) of v n (t) is simply the integral over r.
The three plots (upper row: D 2 , D 3 ; lower left: D 4 ) showing the linear response of v n (t) to the corresponding initial ε n exhibit similar qualitative features.The innermost region of the system -extending up to r R in the case of D 2 , up to r 0.7R for D 3 and D 4 -contributes to v n with the same sign (positive for n = 2 and 4, negative for n = 3) over the whole evolution.Further away from the center come regions that contribute with the opposite sign, not much so for v 2 , more visibly for v 3 and v 4 .A third outer region with the same sign as the innermost one is clearly visible in the case n = 4, and very faintly for n = 3.As time passes by, these regions tend to move towards larger r values, but less markedly than in the similar study for energy-weighted flow [19].
These space-time dependent D n underlie the time dependence of the corresponding v n (t).Thus, the change of sign of v 4 (t) from positive to negative in the upper panels of Fig. 8, with a derivative that turns negative around t/R 1.5, reflects the progressive dominance of the region at intermediate r in the lower left panel of Fig. 10.Similarly, the (small) decrease of v 2 (t) for t/R > ∼ 2 is due to the outer, negatively contributing regions in D 2 (t, r).In the case n = 3 the contributions from the various regions exactly cancel out at every t to yield ∂ t v 3 (t) = 0, while in Ref. [18] an "almost nearly perfect cancellation" resulting in a very small negative v E 3 value was found.Eventually, one can also note that the buildup of the linear D n happens more slowly with increasing n, which possibly reflects the scaling behavior v n (t) ∝ t n+1 in the few collisions regime [32].
The lower right plot of Fig. 10, showing D 4 for the nonlinear response of v 4 to an initial ε 2 , is completely different, with a clear negative contribution at early times and for r < ∼ 1.3R, followed by a positive contribution at later times and for all values of r.In that case one easily checks that the initial eccentricity ε 2 = 0, irrespective of its sign, i.e. the ellipse orientation, generates via the loss term a negative v 4 ∝ −ε 2 2 .Simultaneously, the ellipticity decreases in absolute value, due to the v 2 which is also created.As the negative quadrangular flow v 4 develops, it leads to the development of a positive quadrangularity ε 4 , which is the seed for the positive contributions to v 4 at later times, as seen in the bottom panels of Fig. 8.

IV. SUMMARY
We have investigated anisotropic flow in the fewrescatterings regime in four models, starting with a transport code with elastic binary scatterings, which serves as the reference including all rescatterings in the system.To assess which fraction of the anisotropic flow is carried by particles that escape the system without scattering, we introduced a 2 → 0 version of the code.With the help of a variant of the 2 → 2 code in which particles that have collided once may no longer rescatter, but are accounted for in the final state, we estimate the amount of anisotropic flow at the "single-hit" level.Eventually, we carried out analytical calculations within Boltzmann kinetic theory, including only the loss term of the binary collision kernel and restricting ourselves to linear order in the cross section.Intrinsically the analytical approach and the simulations with the 2 → 0 kernel are unphysical, since energy and momentum are not conserved in the rescatterings.Nevertheless they provide us with a proxy on how much anisotropic flow is created by particles escaping the system without any interaction.
On the other side, the strength of the analytical calculations is that they yield directly a number of known scaling behaviors of the anisotropic flow coefficients, like their dependence on the initial-state eccentricities or their early-time onset, confirming earlier studies [32].Remarkably, the analytical approach at order O(σ) yields v n = 0 for all odd coefficients, but finite values for even ones, which hints at a fundamental difference between odd and even harmonics.Note that we have found elsewhere that odd v n harmonics can be non-zero at order O(σ 2 ) [20].
For even harmonics (v 2 , v 4 , v 6 ), the results of all approaches are very similar when the number of rescatterings in the system is small.In the case of v 4 and v 6 , this holds for both the linear flow response v n ∝ ε n and the nonlinear response like e.g.v 4 ∝ ε 2 2 .The agreement suggests that in the few-rescatterings regime, the even components of the flow signal are to a large extent carried by particles that flew out of the system without colliding, with an anisotropic escape probability reflecting the asymmetric geometry, as advocated for v 2 in AMPT [21].
In contrast, for odd harmonics (v 3 ) the results of the 2 → 2 and 2 → 0 numerical scenarios differ significantly, even in the very few rescatterings regime.Indeed, the former scale roughly linearly with N resc., while the latter rather scale like N 2 resc. .That finding in the 2 → 0 model is consistent with the fact that we find v 3 = 0 in our analytical calculations at order σ -while we found in a parallel study that there is a non-zero v 3 at order σ 2 [32].The v 3 results from the single-hit model also differ significantly from those of 2 → 2 simulations.All in all, the strong dependence of triangular flow on the choice of collision kernel confirms the observation in Ref. [18].In particular, the discrepancy between the approaches demonstrates that in the case of the odd harmonics, the observed v n is not driven by the anisotropic-escape prob-ability, but that the fate of particles after they have undergone a collision does matter.
A clear limitation of the present study is the restriction to a two-dimensional expansion.As we explained in the introduction, this is due to the fact that the small v 3 values require large statistics, which would be too timeconsuming in a three-dimensional study.Indeed, we want to emphasize that previous studies [13,[15][16][17][18][19] of kinetic theory at small opacity relied on solving the (deterministic) Boltzmann equation -with different collision kernels -, while here for the first time10 we used transport simulations at small N resc. .This makes it significantly harder to obtain reliably very small v n values, of order a few 10 −5 at the smallest N resc.we considered (see Figs. 5,8,9).This is even more true in the presence of longitudinal expansion, which dilutes the transverse profile of the system faster, thereby decreasing the anisotropic flow.
That being told, we may still comment on the results one can anticipate in a three-dimensional expansion, in particular a longitudinally boost-invariant one.First, as pointed out in Appendix A 2, the property that v 3 and higher odd harmonics vanish in the loss-term-only calculations is sensitive to the presence of a longitudinal direction: In a three-dimensional geometry, odd harmonics are probably zero at linear order in σ only if the particles are massless and the local momentum distribution in the initial state is independent of position, which is unrealistic.That is, we would anticipate that our finding v 3 ∝ N 2 resc.
may not be robust and be replaced by v n ∝ N resc.for both even and odd n.It is also clear that rescatterings will generally change the longitudinal components of momenta.Thus, it is possible that the agreement we find between all models for even harmonics may not survive the introduction of a third dimension, i.e. that the apparent importance of the anisotropic-escape contribution to the coefficients v 2n may no longer persist.However, we do not see how longitudinal expansion could enhance the effectiveness of the escape mechanism at producing the odd flow harmonics -although it may decrease the relative importance of the component modeled by the gain term of the Boltzmann equation in some regions of phase space.
We would thus conclude that the "escape mechanism" picture cannot account for the whole anisotropic flow signal in systems with very few rescatterings per particle.Within our study, the mechanism is efficient for even harmonics, but not for odd ones.It also means that the details of the (differential) scattering cross section certainly matter for predicting the value of odd anisotropic flow harmonics, as already hinted at by the results on v 3 in Ref. [18], while the even harmonics may be less sensitive.To our knowledge, such a difference in the microscopic "origin" of even and odd flow harmonics has not been reported before in the framework of transport studies. 11ince this will be relevant in systems with small enough multiplicities, our study within a toy transport model (two-dimensional expansion, hard spheres) clearly needs to be replicated with more realistic codes and setups.
where Λ(x) symbolizes the dependence of the local momentum distribution on the position in the transverse plane: this could for instance be a local saturation scale or a local temperature.As previously, we may still assume without loss of generality that G is normalized to unity at every position x when integrating over the whole momentum space -as is e.g. the case if it is a thermal Boltzmann distribution G(p; Λ(x)) ∝ e −|p|/Λ(x) /Λ(x) 2 .The crux is that since G is assumed to be isotropic in momentum space, i.e.only depends on the modulus |p| the normalization of G translates at once into which holds irrespective of whether or not the local momentum distribution is position-dependent.
Let us go back to Eq. ( 7) which gives v n (t) at order O(σ) in the loss-term approach.The free-streaming distributions f f.s. in the integrand are evaluated at positions x − vt resp.x − v 1 t [cf.Eq. ( 5)] that only involve the azimuths of the momenta p, p 1 .The integrals over the moduli |p| and |p 1 | can thus be performed at once using Eq.(A8).That is, effectively the precise form of G -especially its dependence or not on position -does not matter for v n (t) in our analytical approach.Thus, if the odd v n (t) harmonics vanish for a positionindependent initial momentum distribution, as shown in Appendix A 1, then this remains true if G depends on x.
Note that the above proof does not readily generalize to a three-dimensional system nor to massive particles: in such cases, the Møller velocity in the integrand of Eq. ( 7) takes a more complicated form, and in particular it depends on |p| and |p 1 |, so that Eq. (A8) can no longer be used.Similarly, it does not hold either for energyweighted flow coefficients (like v E 3 ), because in that case an extra factor of |p| enters the integrand on the righthand side of Eq. ( 7), which again prevents the use of Eq. (A8).This is consistent with the findings in Ref. [20] in which v 3 is first non-zero at order O(σ 2 ), while v E 3 is already finite at order O(σ).
FIG. 11.Ratio v2(t)/Kn −1 in systems with initially ε2 = 0.15 and on average Nresc.≈ 0.02 (full lines) or 0.14 (dashed) rescatterings per particle for the three scenarios of the transport cascade: 2 → 2 (green), 2 → 0 (blue), single hit (orange).In this Appendix we provide for the sake of reference results for v 2 (Fig. 12) and v 4 (Fig. 13) for systems in which the mean number of rescatterings per particle is about 0.35.For the 2 → 0 resp.single-hit scenarios, this means that approximately 70% of the particles disappear resp.become transparent over the system evolution.Accordingly, the assumption underlying the analytical calculations, that the phase-space distribution f (t, x, p) deviates negligibly at all times from the free-streaming distribution f f.s.(t, x, p) with the same initial condition, is clearly non fulfilled.In addition, it is somewhat clear that if 70% of the particles scatter once in the 2 → 0 or single-hit models, then a significant fraction of them would actually collide several times in the 2 → 2 model: extrapolating the straight-line fit in Fig. 1 indeed gives N resc.≈ 1.67 for a 2 → 2 system with the same initial input Knudsen number as used in the 2 → 0 or single-hit simulations.That is, it is clear from the start that the "full" 2 → 2 and "truncated" 2 → 0 or single-hit systems that lead to N resc.≈ 0.35 are extremely different.
The two plots displaying "linear" flow response, namely v n for an initial non-zero ε n with n = 2 (Fig. 12) or n = 4 (Fig. 13 left), are similar: The results from the simulations with the 2 → 2 collision kernel (green lines) and the 2 → 0 scenario (blue lines) largely differ, by roughly 40% in the case of v 2 , and even yielding signals with opposite signs in the case of v 4 .In contrast, the analytical results are remarkably close to those from the transport calculations with 2 → 2 scatterings, in particular the final values of v 2 or v 4 , which in our view should probably not be over-interpreted.As mentioned in Sect.III C, the agreement for v 4 disappears at higher N resc.values, since the 2 → 2 results become positive.

FIG. 5 .
FIG.5.Time dependence of triangular flow v3 in systems with initially ε3 = 0.15 and on average Nresc.≈ 0.02 (top left), 0.08 (top right) or 0.14 (bottom) rescatterings per particle.The green curves are for systems with elastic binary scatterings, the blue lines for the 2 → 0 collision kernel, and the orange curves for the single-hit model.The constant red line v3 = 0 is the output of the analytical approach.

FIG. 7 .
FIG. 7. Time dependence of energy-weighted triangular flow v E3 in systems with initially ε3 = 0.15 and on average Nresc.≈ 0.02 (top left), 0.08 (top right) or 0.14 (bottom) rescatterings per particle.The green curves are for systems with elastic binary scatterings, the blue lines for the 2 → 0 collision kernel, and the orange curves for the single-hit model.The constant red line v E 3 = 0 is the output of the analytical approach.