Fluctuations of anisotropic flow from the finite number of rescatterings in a two-dimensional massless transport model

We investigate the fluctuations of anisotropic transverse flow due to the finite number of scatterings in a two-dimensional system of massless particles. Using a set of initial geometries from a Monte Carlo Glauber model, we study how flow coefficients fluctuate about their mean value at the corresponding eccentricity, for several values of the scattering cross section. We also show how the distributions of the second and third event planes of anisotropic flow about the corresponding participant plane in the initial geometry evolve as a function of the mean number of scatterings in the system.


I. INTRODUCTION
Anisotropic flow is one of the most prominent observables in heavy-ion collision experiments, especially at high energies. The signal, a modulation in the transverse emission pattern of particles, is believed to be largely due to the lack of cylindrical symmetry of the initial geometry, which is converted by rescatterings into an anisotropy in momentum space [1][2][3].
When the two nuclei collide, the energy they deposit in the region where they overlap is not evenly distributed. The random positions at the instant of the collision of the nucleons, and of the underlying quarks and gluons with their associated color fields, result in an inhomogeneous energy-density profile in the overlap zone. Restricting oneself to the transverse geometry of that zone, described by polar coordinates (r, θ), its azimuthal asymmetries are usually characterized by eccentricities [4][5][6] ε n e inΦn ≡ − r n e inθ r n , where the angular brackets denote an average over the centered energy (or sometimes entropy) density, while the definition given here holds for n ≥ 2. Φ n is the azimuth of the so-called n-th participant plane. In turn, the anisotropies of the transverse momentum distribution of particles are quantified by the successive coefficients of a Fourier series [7] v n e inΨn ≡ e inϕp , with Ψ n the angle of n-th "event plane". Here the brackets denote an average over the momentum distribution of particles in one event, with ϕ p the azimuth of a particle momentum.
As stated above, accounting for the irregularities in the initial geometry of a collision is essential to properly estimate the eccentricity in a given harmonic [4,13]. This randomness in the initial state is often modeled with the help of Monte Carlo (MC) Glauber models [14], which is the approach we use in this paper. In actual measurements of anisotropic flow, the Fourier coefficients are obtained by averaging over many events with varying initial geometries. The event-by-event fluctuations of the eccentricities, which can be studied within a given model of the initial state [15], result in corresponding fluctuations of the flow harmonics within the set of events used for a specific analysis. These flow fluctuations are being extensively investigated, both experimentally [16][17][18][19] and theoretically [10][11][12][20][21][22][23][24][25][26][27].
Most of the phenomenological studies of fluctuations employ hydrodynamics to describe the main stage of the system dynamics, which reflects the overall successful description of bulk observables -especially anisotropic flow -in hydrodynamic studies [1]. However, one could argue that hydrodynamics is possibly too deterministic for a study of fluctuations. Assuming the existence of an underlying particle-based description of the system, hydrodynamics corresponds to the regime in which the (quasi)particles undergo very many scatterings. Thus, hydrodynamics 1 does not capture the possible effect of having only a finite number of collisions per particlewhich could even become quite small for collisions of smaller systems, in which anisotropic flow and its fluctuations are well established [29].
In this article, we focus on the fluctuations of the anisotropic flow harmonics that arise due to the finite number of rescatterings N resc. per particle in a kinetictransport description of the system. We investigate the influence of the average number of particle scatterings on the relationship between elliptic or triangular flow and the corresponding eccentricity. That is, we study how the distribution of v n at a given ε n -i.e. the conditional probability distribution p v|ε (v n |ε n ) -is affected by the average number of rescatterings per particle, extending the exploratory study of Ref. [30].
In the next section we introduce our approach, which consists of two main components: the calculation of an energy density profile from a MC Glauber model and its conversion into a system of massless test particles (Sec. II A) and a two-dimensional transport model for the subsequent evolution of the system (Sec. II B). Section III describes our approach for the analysis of the conditional probability distribution p v|ε (v n |ε n ) as a function of the average number of rescatterings per particle. In Sec. IV we present the results of our statistical analysis in terms of the moments of the conditional probability distribution. We also show how the distributions of event-plane angles relative to the participant-plane angles evolve with the number of rescatterings. Additional results are presented in the Appendices. Eventually we summarize our main findings and discuss them in Sec. V.

II. SETUP OF THE SIMULATIONS
In this section we describe the details of the preparation of the initial states we use, and we introduce the transport code that evolves these initial conditions.

A. Initial state
For the initial states for our transport calculations, we used the TGlauberMC code [31] to generate overlap geometries of two lead nuclei. Using the cross section σ inel NN = (67.6 ± 0.6) mb of inelastic nucleon-nucleon collisions at √ s = 5.02 TeV, we simulated three sets of 10 4 events. The impact parameter b was fixed in each set, at the values b = 0, 6 and 9 fm, respectively, while its direction always lies along the x-axis. From the nucleon positions (x, y) in the transverse plane and the local numbers of binary collisions N coll. (x, y) and of participants N part. (x, y), we first generated an approximate energy density distribution. We used an overall scale factor such that the "hottest spot" with the highest energy density corresponds to a temperature of about 800 MeV. This histogram-like energy density was then smeared as a Gaussian distribution 2 with a width of i.e. approximately 0.7 fm. Note that R N is also the length used in the Glauber model to decide if two nucleons collide. The next step in the generation of our initial condition is the "particlization" of the smooth energy density profile. To that effect, we compute the particle-number density n(x, y) corresponding to the energy distribution, using the formula for a two-dimensional ideal gas of massless particles at thermal equilibrium without fugacity. Particle positions are then sampled from n(x, y) with an acceptance-rejection algorithm. In each event, we sample N p = 5 × 10 5 test particles, irrespective of the actual total energy in the initial state.
To implement local energy conservation as accurately as possible while keeping the running time of the transport code manageable, we employed a coarser grid spacing of 0.1 fm for the particle density grid in the particlization procedure instead of the 0.01 fm used in the energy density calculation from the Glauber model. This allows us to minimize the number of cells on the grid with a sizable energy but no test particle, without having to resort to too large values of N p .
Eventually, we generate the local momentum distribution in the initial state. The energy content of each cell on the grid is evenly distributed between the test particles in the corresponding cell. The particle momenta are rotated by a random angle in the transverse plane to mimic an isotropic initial condition in momentum space. In fact, due to numerical fluctuations even the overall momentum distribution cannot be exactly isotropic in a single event. We ensured the absence of an initial global momentum in each event, i.e. the absence of a directed flow v 1 , but higher-order anisotropies of order N −1/2 p are unavoidable, which we shall discuss in Sec. IV A.
Since we are interested in the dependence on the number of rescatterings in the system of the conditional probability distribution p v|ε (v n |ε n ) describing the anisotropic flow response for a given eccentricity, it would seem preferable to generate initial states with a fixed geometry, instead of our choice. We did not follow this approach for two reasons: On the one hand, we first tried to use "artificial" initial states consisting of distorted Gaussian distributions in the spirit of Ref. [32,33], in which the eccentricities are an input. However we found out that this led to a non-standard v 2 response, namely with a negative cubic term, i.e. K 2,222 < 0 in Eq. (5), while realistic initial conditions yield K 2,222 > 0 [12] (see also  Table VIII). On the other hand, working with a few "realistic" initial conditions, like say a dozen events from the TGlauberMC code, and running the subsequent transport component a large number of times with different seeds for the random-number generator entailed the risk of having picked special initial distributions, that would not reflect the generic behavior. Since the overall dependence on the number of rescatterings that we report appears to be consistent over the whole eccentricity range scanned in fixed-impact parameter MC Glauber events, we think that our approach, albeit not optimal, is valid.

B. Two-dimensional massless transport
The second main ingredient in our simulation is the transport code that performs the time evolution of the initial condition we have just presented. Since the main purpose of our study is to investigate the statistical properties of the physical fluctuations of anisotropic flow arising from the finite number of rescatterings per particle N resc. , we need to simulate a large number of events. To keep the total computational time reasonable, 3 we restricted ourselves to a purely two dimensional system, i.e., our massless test particles live and propagate in the (x, y)-plane only. 4 This obviously means that the absolute flow values that are shown below are not meaningful for phenomenological applications. Nonetheless, since anisotropic flow in heavy ion collisions is to a large extent controlled by transverse dynamics, in particular if one considers flow at midrapidity, we are confident that the qualitative trends with N resc. that we report hereafter are robust and would persist in fully three-dimensional studies.
Our two-dimensional code is based on the covariant transport theory algorithm of C. Gombeaud and J.-Y. Ollitrault described in Ref. [34], to which we refer for further information, in particular on the smooth approach to the hydrodynamic regime. The massless particles we consider are modeled as Lorentz-contracted hard spheres -actually, since the model is two-dimensional and the particles are massless, they appear as (hard) rods with length r -that undergo elastic two-to-two collisions. In the code, the collision times are determined in the laboratory frame and a given scattering is described as a collision between a pointlike particle and a particle of size 2r. In two dimensions, the "hard spheres" modeling implies that total cross section (with dimension of a length) is σ tot = 2r, while an isotropic differential cross section is ensured by relating the scattering angle θ in 3 To be specific, one iteration of the transport algorithm over a single initial condition requires about 213 s CPU time on a 2.2 GHz processor. The parameters chosen for this test are those with the largest number of rescatterings, hence the longest simulation time. 4 That is, in the equilibrated regime the particles obey twodimensional thermodynamics, implying in particular a speed of sound cs = 1/ √ 2.
the center-of-momentum frame of a given collision to the impact parameter d by θ = π(1 − d/r) -which also ensures that the two particles do not rescatter at once after a first collision. The time step of the evolution is chosen such that the particles can collide only once in a time step.
To preserve covariance and locality in the algorithm, and to remain in the regime in which the system dynamics are described by the Boltzmann kinetic theory with a collision kernel including only two-to-two scatterings, it is important that the system should remain dilute. This requirement can be quantified by the dimensionless dilution parameter D defined as the ratio of the mean interparticle distance n −1/2 over the mean free path mfp which has to fulfill D 1. For the actual values of D and of the Knudsen number (9) given in the following, the mean free path is taken as mfp ≡ 1/(nσ tot ), using the average density n of the system in the initial condition of the event under consideration. In all simulations we report hereafter, we found D < ∼ 0.1, which according to Ref. [34] (see in particular Fig. 2) indeed corresponds to the dilute regime, at least as far as anisotropic flow is concerned.
As already stated, we use N p = 5 × 10 5 test particles per event in our simulations. To further decrease the numerical noise on the flow harmonics v n due to the finite number of particles, every Gaussian-smeared initial energy-density profile was converted into a particle system N iter. = 10 times, and the N iter. resulting particle distributions were let to evolve with the transport algorithm independently. Running N iter. events with the same initial energy-density profile and with N p particles each and averaging the flow values over the N iter. iterations is equivalent, from the point of view of the numerical uncertainties on the v n , to considering a single event with N iter. · N p test particles: in either case the resulting uncertainty on a v n is of order (and actually even a factor √ 2 smaller [7]). Performing the N iter. iterations is however faster in practice, since the running time grows faster than linearly with the particle number N p .
For every impact parameter we simulate 10 4 MC energy distributions, each of which gives rise to N iter. initial conditions with N p test particles. In turn, each of these initial particle distributions is used 3 times as input to the transport algorithm, namely with 3 different cross sections σ tot (or equivalently sizes r) for the scatterings of the test particles. As we shall see later, the largest cross section leads to a v 2 or v 3 close to the value in reaches in the hydrodynamic limit, while the other two values correspond to the regime in which v 2 and v 3 vary significantly with σ tot . Instead of using the cross section, we shall characterize the three types of simulations by the Knudsen number Kn, defined as the ratio of the mean free path mfp over a length R characterizing the system size: For the characteristic length R, we use the root mean square radius (with energy-density weighting) of the initial overlap region of the two lead nuclei. That is, R varies event to event, as does the mean free path. Accordingly, each value of σ tot does not correspond to a single Knudsen number, but to a distribution about a mean value Kn . In practice, the relative fluctuation of Kn about Kn is of less than 5 %, so we shall ignore it in the following and only quote the mean value or its inverse. Note that the value of Kn for a given total cross section σ tot does depend on the impact parameter of the lead-lead collision. The Knudsen number, and more precisely its inverse, is related to the average number of rescatterings per particle N resc. in the system. In our simulations, after averaging over events with a fixed impact-parameter value, we find a very simple linear relation as illustrated in Fig. 1. Note that N resc. depends on the system lifetime. In our computations, we let every simulation run for 15 fm/c, irrespective of the impact parameter and thus the characteristic length scale R. We checked that using a longer lifetime does not affect our results for the fluctuations of the anisotropic flow coefficients, although it naturally leads to a larger proportionality coefficient in Eq. (10).  We already mentioned that the three values of the cross section, or equivalently the Knudsen number, that we chose correspond to different regimes for anisotropic flow. This is illustrated in Fig. 2, which shows v 2 and v 3 , divided by the corresponding eccentricity ε 2 or ε 3 , as a function of Kn for simulations at the three impactparameter values. On the same plot we display two curves: First, a fit through the v 2 /ε 2 values with the formula [35] which was shown [34] to provide a good description of a systematical investigation of the behavior of v 2 as a function of Kn. 5 Here v hydro 2 /ε 2 = 0.2621 and A 2 = 1.1967. Secondly, we describe the behavior of v 3 as a function of the number of rescatterings with the higher-order Padé approximation [36] The fit parameters used in Fig. 2 are v hydro 3 /ε 3 = 0.2096, A 3 = 0.6973, B 3 = 1.7922, and C 3 = 13.945.
As one can see, the simulations with the largest number of rescatterings are in the regime where v 2 and v 3 have almost reached their hydrodynamic values. The simulations with the middle Kn value are in the intermediate regime, while those with the largest Kn correspond to the few-collision case. Note that for triangular flow v 3 , whose development necessitates more rescatterings than v 2 [36,37] (although the rise with N resc. seems to be steeper), our simulations with the intermediate value of the cross section are further away from the corresponding hydrodynamic value. In turn, in the simulations with the largest Kn , v 3 is still at the subpercent level, corresponding rather to the far-from-equilibrium regime.

III. QUANTIFYING THE FLOW FLUCTUATIONS
Let us now describe the procedure we adopt to quantify the fluctuations of the anisotropic flow harmonics.
In this theoretical study, we know the initial state of a given event, in particular the orientations Φ n of the symmetry planes of the initial geometry, see Eq. (1), relative to the direction of the impact parameter. Similarly, we can meaningfully compute the flow coefficients and the associated symmetry plane azimuths event by event. More accurately they are computed by averaging over the N iter. iterations with the same initial energy-density profile, which as we argued above yields values of v n reliable at a few 10 −4 level, see Eq. (8).
To exploit this wealth of information, we shall not study the traditional coefficients v n defined in Eq. (2), but rather separately the cosine and sine parts: where · · · denotes an average over particles. Obviously these coefficients obey v n,c = v n cos(nΨ n ) , v n,s = v n sin(nΨ n ), so that we can reconstruct Ψ n from the knowledge of both v n,c and v n,s in an event. Consistent with our use of v n,c and v n,s -which we shall collectively denote v n,c/s -, we characterize the initial state of an event by where · · · now denotes an average weighted with the energy density. These "eccentricities" are clearly related to the real and imaginary parts of ε n defined in Eq. (1). In contrast to the convention for the latter, the coefficients ε n,c/s are not necessarily positive. An advantage of investigating the behaviors of v n,c and v n,s separately is that their geometric causes ε n,c , ε n,s are to a very large extent uncorrelated, as illustrated by the scatter plot of v 2,s (at the end of the evolution) as a function of ε 2,c shown in Fig. 3 for events at b = 6 fm with the largest cross section we simulate. Restricting ourselves to the lowest harmonics n = 2 and 3, we thus already have four independently fluctuating initial-state eccentricities to our disposal, and as many anisotropic flow coefficients in the final state, whose fluctuations we can study as a function of the inverse mean Knudsen number Kn −1 . On the other hand, a drawback is that v n,c and v n,s are not measured in experiments. Yet as we already mentioned, our two-dimensional model is anyway not adequate for quantitative phenomenological studies. 6 At a given impact-parameter value, we simulate 10 4 initial energy density profiles. For each of those, we compute the four eccentricities ε 2,c , ε 2,s , ε 3,c , ε 3,s . Each energy density distribution is converted into N iter. = 10 particle distributions. These are then let to evolve with the transport algorithm, yielding at the end of the evolution flow coefficients v 2,c , v 2,s , v 3,c , v 3,s , which are averaged over the N iter. runs with the same initial condition. Thus, we have to our disposal 10 4 events with known initial eccentricities and final anisotropic flow coefficients, which we can correlate together.
For each cosine or sine harmonic, we remove as a first step events with the smallest and largest eccentricities, 50 each, to get rid of outliers. 7 The remaining events are then divided according to their eccentricity value into five bins with 1950 events for the first and last bins and 2000 events for the central bins. This is illustrated by the different colors of the points in Fig. 4, in which we display v 2,c as a function of ε 2,c for collisions at b = 6 fm for two values of the cross section. We then perform a global fit using either a linear or a linear + cubic dependence, as in Eqs. (3)-(5). The resulting fit function will be denoted byv n,c/s (ε n,c/s ). In practice, the fit with a cubic term was only relevant for the v 2,c and v 2,s at b = 6 and 9 fm. The fit parameters for the four flow coefficients for the various sets of events we consider are listed in Appendix B. We also present in Appendix C the complex initial eccentricity ε 2 e 2iΦ2 and final elliptic flow v 2 e 2iΨ2 for the events of Fig. 4.
A comparison by eye between both panels of Fig. 4 already hints at the larger width of the simulations with more rescatterings (right) compared to those with smaller inverse Knudsen number on the left. To quantify this observation, which also holds for the other flow harmonics, we calculate the moments of the distribution of v n,c/s about the fit function. More precisely, given the values of ε n,c/s and v n,c/s for an event, which we shall denote in the form v n,c/s (ε n,c/s ), we compute powers of the difference v n,c/s (ε n,c/s ) −v n,c/s (ε n,c/s ), where the fit function is evaluated at the same value of the eccentricity. If we had infinite statistics, we would thus be computing the moments of the conditional probability distribution p v|ε (v n,c/s |ε n,c/s ) about the fit function. Since we only have limited statistics, we can only approach the behavior of these conditional probabilities, which is why we use finite-width bins in ε n,c/s . We do not gather all points in a single bin from the start to avoid averaging out opposite trends in different eccentricity intervals-which can indeed happen, as we shall discuss below.
Note that we cannot reasonably estimate p v|ε (v n,c/s |ε n,c/s ) by using the relation in terms of the probability distributions for v n,c/s and ε n,c/s because of the division by very small numbers at the tails of the distributions that it would imply. Instead we tackle the left hand side of Eq. (16) directly. Let us introduce the moments we shall report in Sec. IV, dropping momentarily the c/s indices for brevity.
The first moment about the fit function is where the index i runs over all N events in a given eccentricity bin. We shall not report the value of µ in the following: it is always a small number of order 10 −5 , although it need not vanish in a bin. Yet its average over bins does vanish by construction of the fit function.
The second moment about the fit function, which almost coincides with the variance about the true average in a bin, is defined by A trivial dependence on the inverse mean Knudsen number comes from the fact that the v n values themselves grow (in absolute value) with Kn −1 . To account for this growth, we shall give values of σ 2 v divided by the square K 2 n,n of the linear coefficient from the corresponding fit function.
The third moment, or skewness, is given by Since it involves an odd power it can be either positive or negative. Indeed both possibilities happen as we shall discuss in Sec. IV B, and it is precisely to account for this behavior that all events were not analyzed in a single bin.
Eventually, we shall also present for a fourth-moment related quantity, namely the (excess) kurtosis Note that while the early study [30] only investigated the variance of the conditional probability of elliptic flow v 2 at a fixed eccentricity ε 2 , in this paper we go beyond the Gaussian approximation and study higher moments. We characterize the uncertainty on the value of a given moment M by a variance σ 2 M . The latter are obtained from a delete-d Jackknife algorithm with random sampling [38], where we randomly delete ten percent of the data in a bin and calculate the moments. This is repeated for 3000 random samples and the sample variance is computed.
These moments (18)- (20) -in the case of the skewness, its absolute value |γ 1 | -and their respective variances are finally averaged over the 5 eccentricity bins. For this, we use a weighted average to give less weight to the bins with larger uncertainties. The weighted average M of a given moment M and the associated variance σ 2 M are computed by where the index j runs over the N bins bins, M j is the value of a moment in bin j, and σ 2 M,j the variance σ 2 M in that bin.

IV. RESULTS
Let us now present our results for the behavior of the flow harmonics as a function of the number of rescatterings in the system. In this section we only give a representative sample from our results at the three different impact parameters, b = 0, 6, and 9 fm. Further results are given in the appendices.

A. Initial anisotropic flow
Since the number of particles N p in each event is finite, the initial momentum distribution before any evolution may already show some sizable anisotropic flow: as a rough estimate, one can anticipate that the initial values of v n,c or v n,s should fluctuate about 0 with a dispersion of order 1/ 2N p = 10 −3 [7]. 8 More precisely, these numerical initial-state fluctuations of v n,c/s should even be Gaussian distributed, and the event-plane angle Ψ n should be uncorrelated to the participant-plane orientation Φ n .
To test these expectations, we chose 5000 among the 10 4 Monte Carlo initial conditions that are simulated for each of the three impact parameter values we consider. We computed the corresponding initial flow coefficients v 2,c , v 2,s , v 3,c , v 3,s in each event, which gives us the approximate probability distribution of these harmonicsirrespective of the corresponding spatial eccentricities. As an example, histograms for the probability distributions obtained at b = 0 are shown in Fig. 5.
To characterize these probability distributions, we computed their first moments performing an analysis similar to that described in Sec. III, with a few differences. First, we used a single bin for all events and thus analyze of all them at once irrespective of their eccentricities, retaining all events. Then we directly assumed that the average values of v n,c/s are zero (the actual mean values are of order 10 −5 ) when computing the moments. Finally, in the Jackknife algorithm for the uncertainties on the moments we deleted 500 events from each set, and repeated the calculations over 6000 samples. The resulting variance σ 2 v , skewness (in absolute value) γ 1 , and excess kurtosis γ 2 are given in Table I, in which we actually scaled σ 2 v by its naive estimate 1/2N p for fluctuations of purely numerical origin. We indeed find σ v 1/ 2N p , while the skewness and excess kurtosis are compatible with zero, i.e. the values for Gaussian distributions.  Similar-looking histograms, leading to comparable moments of the probability distributions, are found for b = 6 or 9 fm (see Appendix A), which shows that these initial flow harmonics are driven by numerical fluctuations only.

B. Moments of the final state distribution
We now turn to the flow harmonics in the final state, as a function of Kn −1 . Before we show any result, we emphasize that the trends they exhibit are consistent across different impact parameters. Accordingly, here and in Sec. IV C we only show results at b = 6 fm, postponing those at 0 or 9 fm to Appendices D and E.
In Table II resp. Table III we give the moments (18)- (20) for v 2,c/s resp. v 3,c/s for the three values of the cross section in our simulations. Regarding the variance, here divided by K 2 n,n which monotonously increases with Kn −1 , we observe that it is consistently highest for the  simulations with the largest mean number of rescatterings. On the other hand, the values of σ 2 v in the collisions with the smallest Kn −1 are close to that in the initial state and due to numerical noise (Sec. IV A). The overall trend seems therefore to be that the dispersion of the flow harmonics about their "average" value as represented by the fit functionv n,c/s is increasing with Kn −1 , and that this increase is stronger than that ofv n,c/s itself. This corresponds to the visual impression of a larger width of the scatter plot in the right panel of Fig. 4 mentioned above.
Let us now discuss the skewness. In the tables, we only report absolute values, but γ 1 can be either positive or negative. Before averaging over the eccentricity bins, we actually observed that the sign of γ 1 does depend on the bin, with a clear trend. In bins with a positive ε n,c/s , resulting mostly in a positive final-state v n,c/s , the skewness γ 1 is negative. This signals a distribution with a longer tail towards the left of its positive average, i.e. a distribution skewed towards 0. This behavior was observed for the fluctuations of elliptic flow v 2 in hydrodynamics calculations [26], but as we discuss in next section, the underlying reason probably differs. In turn, in bins with a negative eccentricity we find a positive skewness, which again means a distribution skewed towards 0. To remedy this change of sign of γ 1 , we average the absolute values |γ 1 |. In most cases, the average |γ 1 | decreases with increasing Kn −1 . For the largest cross section we have studied, |γ 1 | seems to be approaching 0, i.e. the value for a Gaussian distribution. On the other hand, in the simulations with the smallest Kn −1 it is far from 0 -and in particular much larger than in the initial state. This hints at the fact that a (very) small number of rescatter-ings 9 would lead to a v n,c/s distribution with a non-zero mean and a strong skewness towards zero, and that |γ 1 | then decreases with increasing Kn −1 .
Regarding the excess kurtosis γ 2 , the rather large error bars make it difficult to draw definite conclusions. If one nevertheless looks for overall trends, one can identify two. First, γ 2 seems to favor positive values, and often values that are not compatible with 0. If we try to be more specific, one could say that γ 2 is always positive for v 2 , while it is at times negative (yet compatible with 0 with error bars) for v 3 , which happens to coincide with the findings in Ref. [27] in a hydrodynamic approach. A second possible trend is a non-monotonic behavior with Kn −1 , with a positive maximum -signaling a distribution that is more peaked than a Gaussian -at the intermediate value that we simulated, while for the smallest and largest Kn −1 the behavior is more Gaussian.
In general we can say that for many collisions the absolute value of the skewness and the kurtosis approach the values for an almost Gaussian shape, as in the initial state. Regarding the variance of the final-state flow harmonics, it is larger than in the initial state, Eventually, the general effect of considering a finite number of rescatterings seems be a larger absolute skewness than in the hydrodynamic case, as well as a larger kurtosis.

C. Event-plane angle distributions
From the cosine and sine parts of a flow harmonic, we can calculate the corresponding flow angle Ψ n . Similarly from ε n,c and ε n,s we extract the orientation Φ n of the nth symmetry plane in the initial state. Therefore, we can study probability distributions p(n(Ψ n − Φ n )) for n = 2 and 3. These are displayed in Fig. 6 for our simulations at b = 6 fm, as a function of Kn −1 . The results for the simulations at 0 or 9 fm are in Figs. 11, 12 in Appendix E. In the plots, we also show the distributions about the respective participant plane of the flow angles Ψ 2 and Ψ 3 in the initial state of our simulations, with the label Kn −1 = 0. These initial distributions are flat, up to fluctuations, which shows that there is no correlation between the initial flow, due to numerical noise, in our simulations and the geometry of the events.
Regarding the distributions of Ψ n − Φ n in the case when the flow angle Ψ n is that at the end of the evolution, let us discuss first the results for n = 3. At the smallest value of Kn −1 ≈ 0.3-0.4, a still rather broad distribution develops with a maximum at Ψ 3 − Φ 3 = 0. This clearly signals the onset of the transfer from the initial anisotropy in position space into a signal in momentum space. At the 10 times larger Kn −1 , the distribution is now very peaked around Ψ 3 − Φ 3 = 0: the flow angle is almost aligned along the symmetry plane of 9 Unfortunately, smaller than in the simulations we made. the initial geometry. Eventually, the distribution at the largest value of Kn −1 is less peaked, with thicker tails at larger values of |Ψ 3 − Φ 3 |. We shall come back in next section to this behavior, which was not anticipated but is consistently observed at all three values of the impact parameter, and also in the second harmonic as we discuss now.
Coming to the case n = 2, let us recall our comment at the end of Sec. II B regarding the slower onset of v 3 compared to v 2 as a function of Kn −1 [36,37]. This behavior, which was known for the values of the flow coefficients, is also visible at the level of the flow angles. Indeed, at the smallest value of Kn −1 , we find that the distribution of Ψ 2 − Φ 2 is already very peaked around 0, much more so than p(3(Ψ 3 −Φ 3 )). This is true not only at b = 6 or 9 fm, in which the initial-state geometry yields a clear marked Φ 2 , but even (although to a lesser extent) at b = 0 fm. As the number of rescatterings in the system or equivalently Kn −1 increases, the peak persists. Yet at the highest Kn −1 , we again observe some broadening of the distribution, as in the case n = 3.
For the values of Kn −1 for which the peak is clearly marked, its width is of the same magnitude as reported in hydrodynamic [21] or hybrid [20] simulations.
Finally, let us mention that none of the distribution p(n(Ψ n − Φ n )) can be reasonably fitted by a (truncated) Gaussian nor by a beta distribution.

V. DISCUSSION
We have investigated the fluctuations of anisotropic flow in a two-dimensional system of massless particles undergoing elastic binary scatterings, as a function of the mean number or collisions per particle -which we vary by choosing various values of the scattering cross section. More specifically, we study how each flow coefficient fluctuates about the "preferred" value for the corresponding eccentricity, for instance, how the cosine coefficient v 2,c fluctuates event-to-event about the valuē v 2,c (ε 2,c ) defined by the best fit to the scatter plot of v 2,c vs. ε 2,c with Eq. (5) or its linear version. These fluctuations are characterized by successive moments (variance, skewness, excess kurtosis). In addition, we studied the distribution of the event plane Ψ 2 or Ψ 3 about the corresponding participant plane Φ 2 resp. Φ 3 in the initial state.
Let us gather the salient findings that appear to be robust across the different impact-parameter values we have studied.
• The variance σ 2 v of the fluctuations of v n,c/s about the mean valuev n,c/s (ε n,c/s ), scaled by the overall size of the flow harmonic, first decreases with Kn −1 and increases again.
• The distributions of v n,c/s are always skewed towards 0. The skewness first jumps to a sizable value as flow sets on -to fix ideas, |γ 1 | is of order 0.5 -, then decreases again as the number of rescatterings increases.
• The excess kurtosis γ 2 is predominantly positive, i.e. the distribution is more peaked than a Gaussian with the same variance. It tends towards 0 with increasing Kn −1 .
• The distribution of event plane angles about the orientation of the corresponding participant plane are peaked. The peak broadens for the largest value of Kn −1 .
It seems to us that a consistent picture emerges from these features. First, let us emphasize that the fluctuations we discussed above are not the usual fluctuations arising from those of the eccentricity in a centrality bin. They are fluctuations at fixed eccentricity -which is why we referred to the conditional probability distribution p v|ε (v n,c/s |ε n,c/s ) -, whose moments are averaged over the range of eccentricities found in collisions at a given impact parameter, to increase their significance. That is, we are looking at fluctuations in the response function of anisotropic flow to the initial eccentricity [15]. The physical origin of these fluctuations is essentially the finite number of scatterings that build up the anisotropic flow from the initial eccentricity of the geometry. In turn, the value about which a given flow coefficient fluctuates is that which would be obtained when solving the (fully deterministic) kinetic Boltzmann equation at the same Knudsen number or equivalently opacity [33,37]. 10 Similarly, the fluctuations we have studied are absent from 10 Here we implicitly assume that the Boltzmann equation in terms of a smooth single-particle phase-space distribution is the actual evolution equation describing the system dynamics, which should hold in the dilute regime.
hydrodynamic simulations, at least of those implementing non-fluctuating hydrodynamics.
That being told, how can we interpret the trends in the moments of about the mean value? The first feature is that the fluctuations from the finite number of rescatterings are non-Gaussian, meaning that the central limit theorem does not apply, precisely because there are not enough scatterings. The generic skewness towards zero seems to be intuitively clear, and comes from the fact that the initial value of each flow coefficient is zero and that v n,c/s is bounded in absolute value. On the other hand, the origin of the peaking of the distribution encoded in the kurtosis is not clear to us. Yet that this positive excess kurtosis tends to disappear in the collisions with the largest cross section seems to us to be correlated to the second main feature we observed, namely the broadening of the Ψ n − Φ n distributions in the same regime, which we also find consistently but do not really understand. A possibility -that we have not further tested -is that for the largest Kn −1 nonlinear "contamination" from other harmonics start to be sizable, as e.g. a contribution ε 2 ε 4 to (the fluctuations of) v 2 , both in value and orientation.
The present study can in natural way be extended in several directions. First, one can think of looking at simulations with either less or more rescatterings per particle. Less rescatterings, i.e. a larger Knudsen number, would allow one to study the onset of the skewness, or of the alignment of the event plane angle Φ 2 along the participant plane orientation Ψ 2 . The issue is that less rescatterings mean a smaller flow signal, so one needs to increase the statistics to recover meaningful results. In the other limit, simulations with a smaller Knudsen number would allow to probe the "hydrodynamic limit" more carefully: do the fluctuations of v n,c/s about their mean value due to the finite number of scatterings become smaller -which would be our expectation? How do the distributions of Ψ n − Φ n evolve? The issue with those simulations is that one must ensure that the system remains in the dilute regime, so that cranking up the cross section is not enough, but should be accompanied by a decrease in the number of test particles. In turn, this has to be compensated by an increase in the number of simulated events, to recover small error bars.
Another direction is to look at higher flow harmonics like v 4 or v 5 , or at the correlations between the fluctuations of different harmonics [17]. Here again the signals will be smaller, thus statistics-costly. In addition, these higher harmonics are significantly affected by lower-order eccentricities [6,10,32,39,40], which also begs for more statistics to disentangle the various influences. Eventually, instead of looking only at the average flow coefficients, one can investigate how transverse momentum dependent v n (p T ) are affected by the number of rescatterings.
Going beyond our fixed impact-parameter study, one needs to turn to fixed-centrality studies, in which the fluctuations due to the finite number of rescatterings combine with eccentricity fluctuations. This is probably the simplest extension of our work, and a first step towards investigations that can meaningfully be used for phenomenological comparisons with experimental data. Anticipating on these more realistic studies, one can reasonably imagine that the convolution of the dispersion of the eccentricities within a centrality bin on the one hand, and the fluctuation of the anisotropic-flow response due to the final number of rescatterings will lead to a distortion of the non-Gaussianities of the flow-fluctuation pattern compared to the underlying eccentricity fluctuations. In particular, since we consistently found that the distribution of a flow harmonic at fixed eccentricity is skewed towards zero, one can expect that the distributions of v 2 or v 3 in a centrality bin should be more skewed towards zero than the distributions of ε 2 or ε 3 . Since the average number of rescatterings decreases as the collisions become more peripheral, this "extra negative skewness" of the probability distribution of v n should be more pronounced in the peripheral bins than in the central ones, although a quantitative statement necessitates more realistic computations.   Fig. 7 for the anisotropic flow harmonics in the initial state computed for 5000 events at b = 6 fm. The variance σ 2 v is divided by the rough estimate 1/2Np.    In Tables VI, VII, and VIII we gather the parameters of the fit functions to the scatter plots of v n,c/s (in the final state) as function of ε n,c/s . Linear fits are used throughout, except for v 2,c/s at b = 6 and 9 fm, where the nonlinear ansatz (5) yields better fits.
The fit parameters are consistently very similar for cosine and sine harmonics under the same conditions of b and Kn −1 . All fit parameters increase with Kn −1 at a given impact parameter. This is of course long known for the linear parts of the fits, but also holds for the nonlinear parameters K 2,222 , as recently reported in Ref. [37].  In this Appendix, we show for 10 4 events at b = 6 fm the dispersion of the complex eccentricity ε 2 e 2iΦ2 in the initial state ( Fig. 9 -remember that the impact parameter lies at Φ 2 = 0) and of the final state elliptic flow v 2 e 2iΨ2 for the smallest and largest values of the inverse Knudsen number (Fig. 10). We also indicate the five bins in ε 2,c for those events (see also Fig. 4). In particular, one sees that they do not map onto similar bins in v 2,c = Re v 2 e 2iΨ2 , due to the fluctuation of the response. Figure 9 also shows that the events with "outlying" ε 2,c that we remove from our statistical analysis of v 2,c do not have at the same time an "outlying" ε 2,s , and thus are not the same as are left aside in the study of v 2,s . In this Appendix, we list the moments of the conditional probability distributions p v|ε (v n,c/s |ε n,c/s ) about the mean valuev n,c/s (ε n,c/s ) for our sets of events at b = 0 (Tables IX and X) and 9 fm (Tables XI and XII)  In this Appendix, we show the distributions of each event-plane angle Ψ n (in the final state) about the corresponding participant-plane orientation Φ n (in the initial state) for our sets of events at b = 0 (Fig. 11) and 9 fm (Fig. 12). --π/2 --π/4 0 π/4 π/2 2(Ψ 2 --Φ 2 ) --π/2 --π/4 0 π/4 π/2 2(Ψ 2 --Φ 2 )