# Even anisotropic-flow harmonics are from Venus, odd ones are from Mars

Benedikt Bachmann  Nicolas Borghini<EMAIL_ADDRESS>Nina
Feld<EMAIL_ADDRESS>Hendrik Roch<EMAIL_ADDRESS>Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany

###### Abstract

We test the ability of the “escape mechanism” to create the anisotropic flow
observed in high-energy nuclear collisions. We compare the flow harmonics
$v_{n}$ in the few-rescatterings regime from two types of transport
simulations, with $2\to 2$ and $2\to 0$ collision kernels respectively, and
from analytical calculations neglecting the gain term of the Boltzmann
equation. We find that the even flow harmonics are similar in the three
approaches, while the odd harmonics differ significantly.

##  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 [7, 8, 1, 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 part of the measured anisotropic flow signal, namely the
elliptic flow $v_{2}$, could to a large extent arise not from the many
rescatterings of the produced particles, but 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/\psi$ quarkonia [24].

In this paper, we want to further explore the anisotropic escape scenario, and
in particular study whether its predictions for the flow harmonics are similar
to those within a more classical picture where collisions are essential for
creating anisotropic flow. We focus on the few-rescatterings regime, in which
we can extend our previous analytical calculations which we complement with
numerical transport simulations (Sect. II). 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.

##  II Methods

To investigate the importance of the “escape mechanism” for the production of
anisotropic flow when particles undergo very few rescatterings, we perform
three types of calculations. On the one hand, numerical simulations with a
transport code, with two different collision kernels: first a $2\to 2$ kernel
that implements elastic binary collisions, and gives as reference the “total”
anisotropic flow produced in a semi-realistic system. Secondly, a $2\to 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\to 0$
scenario.

We begin with introducing the analytical approach (Sect. II.1), together with
the initial conditions we use for both analytical and numerical calculations.
We then briefly present our transport setups (Sect. II.2). 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 $\hbar=c=1$, while
$(r,\theta)$ are polar coordinates in the transverse plane, with their origin
at the center of the system in its initial state.

###  II.1 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

$p^{\mu}\partial_{\mu}f(t,{\bf x},{\bf p})={\cal C}[f(t,{\bf x},{\bf p})].$
(1)

Instead of a full collision kernel with detailed balance, we shall only
consider the loss term of binary scatterings

${\cal C}_{\rm loss}[f(t,{\bf x},{\bf p})]=-\frac{E_{\bf
p}}{2}\\!\int\\!f(t,{\bf x},{\bf p})f(t,{\bf x},{\bf p}_{1})v_{\rm
rel.}\sigma\,{\rm d}^{2}{\bf p}_{1},$ (2)

with $E_{\bf p}$ the energy of the particle with momentum ${\bf p}$, $v_{\rm
rel.}$ the Møller velocity, and $\sigma$ the total cross section. Note that
this implies that particle number is not conserved in the evolution. For
massless particles in two dimensions, $v_{\rm rel.}=1-\cos(\varphi_{\bf
p}-\varphi_{1})$ where $\varphi_{\bf p}$ resp. $\varphi_{1}$ is the azimuthal
angle of momentum ${\bf p}$ resp. ${\bf 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

$v_{n}(t)=\frac{\displaystyle\int\\!f(t,{\bf x},{\bf p})\cos(n\varphi_{\bf
p})\,{\rm d}^{2}{\bf x}\,{\rm d}^{2}{\bf p}}{\displaystyle\int\\!f(t,{\bf
x},{\bf p})\,{\rm d}^{2}{\bf x}\,{\rm d}^{2}{\bf p}},$ (3)

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:

$\displaystyle\partial_{t}v_{n}(t)=$ $\displaystyle\
\frac{1}{N(t)}\\!\int\\!\partial_{t}f(t,{\bf x},{\bf p})\cos(n\varphi_{\bf
p})\,{\rm d}^{2}{\bf x}\,{\rm d}^{2}{\bf p}$ (4)
$\displaystyle-\frac{\partial_{t}N(t)}{N(t)}v_{n}(t).$ (5)

Using the Boltzmann equation to replace $\partial_{t}f$ in the integrand, the
term involving the spatial gradient of $f$ gives zero after integrating over
${\bf x}$, since $f$ vanishes at infinity. There remains only the contribution
from the collision term, which at leading order is a priori linear in
$\sigma$. 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. (5),
i.e. approximate $N(t)\simeq 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,{\bf x},{\bf p})$ by the free-streaming
distribution $f_{\rm f.s.}(t,{\bf x},{\bf p})$ that coincides with $f$ in the
initial state [22, 23, 13]:

$f_{\rm f.s.}(t,{\bf x},{\bf p})=f^{(0)\\!}({\bf x}-{\bf v}t,{\bf p}),$ (6)

where ${\bf v}\equiv{\bf p}/|{\bf p}|$ while $f^{(0)}({\bf x},{\bf p})$
denotes the initial distribution (at $t=0$), to which we come back hereafter.

In the second line of Eq. (5), $\partial_{t}N(t)$ is of order ${\cal
O}(\sigma)$ (or higher). In absence of initial anisotropic flow in the system,
$v_{n}(t)$ is also of order ${\cal O}(\sigma)$, so that the whole term is at
least quadratic in $\sigma$: accordingly, we shall neglect it hereafter. Note
however that this term contributes at linear order in $\sigma$, and thus may
not be dropped, if there is some anisotropic flow in the initial state.

All in all, we replace the evolution equation (5) with

$\partial_{t}v_{n}(t)=\frac{1}{N}\\!\int\\!\frac{{\cal C}[f_{\rm f.s.}(t,{\bf
x},{\bf p})]}{E_{\bf p}}\cos(n\varphi_{\bf p})\,{\rm d}^{2}{\bf x}\,{\rm
d}^{2}{\bf p},$ (7)

valid to linear order in $\sigma$, irrespective of the choice of collision
term — as long as the latter is ${\cal O}(\sigma)$. Inserting the loss term
(2) as collision kernel and integrating over time yields

$\displaystyle v_{n}(t)$
$\displaystyle=-\frac{\sigma}{2N}\\!\int_{0}^{t}\\!\int\\!f_{\rm
f.s.}(t^{\prime},{\bf x},{\bf p})f_{\rm f.s.}(t^{\prime},{\bf x},{\bf
p}_{1})[1-\cos(\varphi_{\bf p}-\varphi_{1})]\cos(n\varphi_{\bf p})\,{\rm
d}^{2}{\bf x}\,{\rm d}^{2}{\bf p}\,{\rm d}^{2}{\bf p}_{1}\,{\rm
d}t^{\prime}+{\cal O}(\sigma^{2})$ (8)
$\displaystyle\equiv\int_{0}^{t}\\!\int_{0}^{\infty}\\!D_{n}(t^{\prime},r)\,r\,{\rm
d}r\,{\rm d}t^{\prime}+{\cal O}(\sigma^{2}),$ (9)

where the second line defines the angle-averaged local production rate of
$v_{n}$ [18], which we shall discuss in Sect. III.5. Note that in these
expressions we explicitly assumed $v_{n}(t=0)=0$ in the initial state.

In our analytical approach the flow coefficients (9) depend directly on the
initial phase space distribution $f^{(0)}$ via Eq. (6). 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:

$f^{(0)}({\bf x},{\bf p})=F({\bf x})G({\bf p}),$ (10)

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 (11). 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({\bf p})$ [25], but lead to lengthier expressions for the flow
coefficients — whose evolution at linear order in $\sigma$ is no longer
governed by Eq. (7) as mentioned above.

In position space, we choose as initial density a distorted Gaussian
distribution111In Appendix C we briefly present results using an alternative
initial density.

$F(r,\theta)=\frac{N{\rm e}^{-r^{2\\!}/2R^{2}}}{2\pi
R^{2}}\Bigg{[}1-\sum_{j=2}^{6}\tilde{\varepsilon}_{j}{\rm
e}^{-r^{2\\!}/2R^{2}}\bigg{(}\frac{r}{R}\bigg{)}^{\\!\\!j}\cos(j\theta)\Bigg{]},$
(11)

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 [15, 14, 26, 16, 17, 18, 19], as it allows
one to introduce at will in the initial state different and independent types
of “eccentricities” [27, 28, 29]

$\displaystyle\varepsilon_{n}{\rm e}^{{\rm i}n\Phi_{n}}\equiv-\frac{\langle
r^{n}{\rm e}^{{\rm i}n\theta}\rangle}{\langle r^{n}\rangle},$ (12)

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 (11) yields at once $\Phi_{n}=0$ — which we may assume without loss
of generality since we shall always consider only a single non-zero
$\varepsilon_{n}$ at a time — and

$\varepsilon_{n}=\frac{(n-1)!}{2^{\frac{2+n}{2}}\Gamma(\frac{n}{2})}\tilde{\varepsilon}_{n},$
(13)

that is for the first harmonics $\varepsilon_{2}=\tilde{\varepsilon}_{2}/4$,
$\varepsilon_{3}=\tilde{\varepsilon}_{3}/\sqrt{2\pi}$,
$\varepsilon_{4}=3\tilde{\varepsilon}_{4}/4$, and so on. Note that the
parameters $\tilde{\varepsilon}_{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, $\tilde{\varepsilon}_{n}$ should be such
that $\varepsilon_{n}$ remains smaller than $\varepsilon_{n,\text{max}}\simeq
0.35$. In our calculations, both analytical and numerical, we choose
$\tilde{\varepsilon}_{n}$ such that $\varepsilon_{n}=0.15$ or smaller.

###  II.2 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_{\rm p}$ Lorentz-
contracted hard spheres — or rather hard disks, since they are two-dimensional
— with radius $(N/N_{\rm p})\sigma/2$, where $\sigma$ 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_{\rm p}$ and $\sigma$ 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 $\ell_{\rm mfp}$.

For the simulations with the $2\to 0$ collision kernel we use the same
transport algorithm as in the $2\to 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\to 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.

In Ref. [21] the authors used a similar approach with $2\to 2$ and $2\to 0$
collisions. The difference to our $2\to 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].

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 (11), 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_{\rm p}$ ranging between $2\times
10^{5}$ and $2\times 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_{\rm 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 $\sqrt{N_{\rm iter.}}$. Since the simulation time
grows with $N_{\rm p}^{3/2}$, performing multiple iterations with less test
particles is computationally less costly than performing a single simulation
with $N_{\rm iter.}N_{\rm p}$ particles.222In our simulations, $N_{\rm
iter.}N_{\rm p}$ is always larger than $10^{9}$.

Starting from Eq. (11), the average particle-number density per unit surface
is $N/4\pi R^{2}$. Using the latter to define a mean free path $\ell_{\rm
mfp}$, we quantify the rarity or abundance of rescatterings by the Knudsen
number

${\rm Kn}\equiv\frac{\ell_{\rm mfp}}{R}=\frac{4\pi R}{N\sigma}.$ (14)

Indeed, the average number of rescatterings per particle $N_{\rm resc.}$ over
the whole evolution — in practice, until $t/R=30$ in the simulations — nicely
scales with ${\rm Kn}^{-1}$ in the $2\to 2$ scenario, see Fig. 1.

Figure 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. (14), for the $2\to 2$ (green squares, fit with $N_{\rm resc.}\approx
0.529\,{\rm Kn}^{-1}$) and $2\to 0$ (blue circles, fitted with a quadratic
ansatz: dashed line) scenarios. The red line $N_{\rm resc.}={\rm
Kn}^{-1}/\sqrt{\pi}$ is the prediction of the analytical approach.

In the $2\to 0$ model, a given $N_{\rm resc.}$ requires a larger ${\rm
Kn}^{-1}$, i.e. cross section, since particles disappear over the evolution
and thus can never scatter twice. In the following we mostly present results
for $N_{\rm resc.}\approx 0.02$, well in the few-rescatterings regime, and
0.14 — for which the approximation becomes less justified —, as well as
$N_{\rm resc.}\approx 0.35$ in Appendix B.

##  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 (11).
Calculations with a slightly different initial profile, whose results are in
qualitative agreement with the findings of this Section, are given in Appendix
C.

###  III.1 Elliptic flow

Figure 2: Top: Time dependence of elliptic flow $v_{2}$ in systems with
initially $\varepsilon_{2}=0.15$ and on average $N_{\rm resc.}\approx 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\to 0$
scenario and the red ones show the analytical result (16). The bottom panels
show ratios of the curves from the upper panels.

Let us start with elliptic flow $v_{2}$ [6]. As initial geometry we consider
the profile (11) with all $\tilde{\varepsilon}_{j}=0$ except for
$\tilde{\varepsilon}_{2}$, chosen such that $\varepsilon_{2}=0.15$ (up to
numerical fluctuations in the simulations). The time dependence of $v_{2}$ in
the transport approach is shown in Fig. 2 with both $2\to 2$ (green) and $2\to
0$ (blue) collision kernels, for $N_{\rm resc}\approx 0.02$ (left panel) and
$N_{\rm resc}\approx 0.14$ (right panel).333Results in a system with $N_{\rm
resc.}\approx 0.35$ are shown in Fig. 8. At $t=0$ we indicate as an error bar
the typical value $1/\sqrt{2N_{\rm iter.}N_{\rm p}}$ of $v_{2}$ induced by
numerical fluctuations in the initial state. We also show in red the result
from the analytical calculation, namely

$\displaystyle v_{2}(t)=\frac{8}{27}{\rm Kn}^{-1\,}\varepsilon_{2\,}{\rm
e}^{-2t^{2}/3R^{2}}$
$\displaystyle\bigg{[}\\!\bigg{(}\\!\frac{3R}{t}+\frac{2t}{R}\\!\bigg{)}I_{1}\bigg{(}\\!\frac{2t^{2}}{3R^{2}}\\!\bigg{)}$
(15)
$\displaystyle\quad-\frac{t}{R}I_{0}\bigg{(}\\!\frac{2t^{2}}{3R^{2}}\\!\bigg{)}\bigg{]}$
(16)

with modified Bessel functions of the first kind $I_{0}$ and $I_{1}$, and
where the value of ${\rm Kn}^{-1}$ is chosen such that it gives the same
$N_{\rm resc.}$ as in the numerical calculations. Note that Eq. (16) yields
$v_{2}(t)\propto t^{3}$ at early times $t\ll R$, as pointed out in previous
studies [12, 31, 23, 32].

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

$v_{2}(t)\sim\frac{\sum_{k=3}^{5}a_{k}(t/R)^{k}}{1+\sum_{k=1}^{5}b_{k}(t/R)^{k}}$
(17)

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\gtrsim 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\to 0$ scenario either to those of the $2\to 2$ case or to the
analytical value (16), and show these ratios in the narrow lower panels in
Fig. 2.

The profiles of $v_{2}(t)$ are similar in the three approaches, with a slow
onset, followed by an almost linear rise, that eventually saturates. $v_{2}$
reaches its maximum value for $t/R\approx 2$, and decreases a little
afterwards, barely in the $2\to 2$ scenario. 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\to 0$ model differs at most by 20% from that in the “full” $2\to 2$
case. In addition, the results of the $2\to 0$ scenario are extremely well
reproduced by the analytical calculations accounting only for the loss term
for $N_{\rm resc.}\approx 0.02$. The agreement is less impressive but still
very good at the larger $N_{\rm resc.}$, which is easily understood: The
analytical results are derived at linear order in the cross section, or
equivalently $\mathrm{Kn}^{-1}$. As the latter increases, higher order
contributions to $v_{2}$, which are always present in the $2\to 0$
simulations, become more sizable, and lead to the departure between the
analytical results and the $2\to 0$ values. Indeed we have shown in Ref. [20]
— yet only for early times — that pushing the analytical calculation to higher
order in $\sigma$ improves the agreement with the $2\to 0$ results. In
contrast, for $t/R\geq 2$, when fewer collisions take place, the results of
both approaches are again very parallel.

Thus, 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].

###  III.2 Triangular flow

We turn next to triangular flow $v_{3}$ [27], using now an initial geometrical
profile (11) with only a non-zero $\tilde{\varepsilon}_{3}$, such that
$\varepsilon_{3}=0.15$. The results of our various calculations for the time
dependence of $v_{3}$ are displayed in Fig. 3, for systems with $N_{\rm
resc.}\approx 0.02$ (top left) 0.08 (top right) or $0.14$ (bottom).

Figure 3: Time dependence of triangular flow $v_{3}$ in systems with initially
$\varepsilon_{3}=0.15$ and on average $N_{\rm resc.}\approx 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\to 0$ collision kernel. The constant red line $v_{3}=0$ is the output of
the analytical approach.

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.5, 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
$\sigma$, 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.444The analytical
calculation with an initial momentum anisotropy leads to a non-vanishing and
slightly evolving $v_{3}$. The latter is however negligible compared to the
values of the numerical simulations and therefore not shown in Fig. 3.
Secondly, the simulations within the $2\to 0$ model give a clear non-zero
signal, in contrast to the analytical result. This hints that in the $2\to 0$
simulations, which include all orders in the cross section, $v_{3}$ arises at
a higher order in $\sigma$.

Thirdly, the results of the $2\to 0$ scenario clearly do not resemble those of
the $2\to 2$ model. For $N_{\rm resc.}\approx 0.02$, the $2\to 0$ results lie
about a factor 1.5 below, while they are larger for $N_{\rm resc.}\approx
0.08$ and 0.14. Actually, $v_{3}$ scales approximately with $N_{\rm resc.}$ in
the $2\to 2$ simulations, but rather with $N_{\rm resc.}^{2}$ in the $2\to 0$
scenario.

Figure 4: Time dependence of energy-weighted triangular flow $v_{3}^{E}$ in
systems with initially $\varepsilon_{3}=0.15$ and on average $N_{\rm
resc.}\approx 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\to 0$ collision kernel. The
constant red line $v_{3}^{E}=0$ is the output of the analytical approach.

Several recent studies investigated the “energy weighted triangular flow”
$v_{3}^{E}$, 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. 4, $v_{3}^{E}$ — computed in the same systems as used
for Fig. 3 — again differs a lot in the $2\to 2$ and $2\to 0$ scenarios. This
is especially true at times $t\gtrsim R$. At earlier times, the results with
the two collision kernels 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.1, we find
that for $v_{3}$ the $2\to 0$ scenario differs significantly from the $2\to 2$
model. In parallel, the triangular flow from the analytical approach
considering only the loss term at first order in $\sigma$ is also at variance
with the results of numerical simulations. This 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.

###  III.3 Quadrangular flow

With quadrangular flow $v_{4}$, the situation is again simpler than for
$v_{3}$. Anticipating on what shall now present, the overall trend is the same
as for elliptic flow $v_{2}$: the results of the numerical $2\to 2$ and $2\to
0$ 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.

Figure 5: Time dependence of quadrangular flow $v_{4}$ in systems with on
average $N_{\rm resc.}\approx 0.02$ (left) or $0.14$ (right) rescatterings per
particle. The top panels are for a system with an initial quadrangularity
$\varepsilon_{4}=0.15$ and $\varepsilon_{2}=0$; the bottom panels for an
initial state with $\varepsilon_{2}=0.15$ and $\varepsilon_{4}=0$. The green
curves are for systems with elastic binary scatterings, the blue lines for the
$2\to 0$ scenario and the red ones show the analytical results (18) (top) or
Eq. (19) (bottom).

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 $\varepsilon_{n}$, but also due to nonlinear effects mixing other
eccentricities.555To be more accurate, according to our present knowledge the
lower flow harmonics $v_{n}$ with $n\leq 3$ are only minimally affected by
such nonlinear effects involving eccentricities $\varepsilon_{k}$ with $k\neq
n$. Thus, $v_{4}$ may be caused not only by the “quadrangularity”
$\varepsilon_{4}$, but also by the “ellipticity” $\varepsilon_{2}$ [33, 29,
34, 35, 14].

Indeed, our analytical calculation for $v_{4}$ assuming only a non-vanishing
$\tilde{\varepsilon}_{4}$ in the initial state yields

$v_{4}(t)=\frac{16}{1215}{\rm Kn}^{-1\,}\varepsilon_{4\,}{\rm
e}^{-2t^{2}/3R^{2}}\bigg{[}\bigg{(}\frac{162R^{3}}{t^{3}}+\frac{63R}{t}+\frac{24t}{R}+\frac{5t^{3}}{R^{3}}\bigg{)}I_{1}\bigg{(}\frac{2t^{2}}{3R^{2}}\bigg{)}-\bigg{(}\frac{54R}{t}+\frac{21t}{R}+\frac{5t^{3}}{R^{3}}\bigg{)}I_{0}\bigg{(}\frac{2t^{2}}{3R^{2}}\bigg{)}\bigg{]}.$
(18)

Assuming instead that only a non-vanishing $\tilde{\varepsilon}_{2}$ is
initially present, we obtain

$\displaystyle v_{4}(t)=-\frac{1}{10}{\rm Kn}^{-1\,}\varepsilon_{2\,}^{2}{\rm
e}^{-t^{2}/R^{2}}\bigg{[}\bigg{(}\frac{48R^{3}}{t^{3}}+\frac{28R}{t}+\frac{16t}{R}+\frac{5t^{3}}{R^{3}}\bigg{)}I_{1}\bigg{(}\frac{t^{2}}{R^{2}}\bigg{)}-\bigg{(}\frac{24R}{t}+\frac{14t}{R}+\frac{5t^{3}}{R^{3}}\bigg{)}I_{0}\bigg{(}\frac{t^{2}}{R^{2}}\bigg{)}\bigg{]}.$
(19)

Obviously, the terms on the right hand sides of these equations add up if the
initial state contains both $\tilde{\varepsilon}_{2}$ and
$\tilde{\varepsilon}_{4}$. These analytical results are compared to those of
numerical simulations with both $2\to 2$ (green) and $2\to 0$ (blue) collision
kernels in Fig. 5: the plots in the top panels are with
$\tilde{\varepsilon}_{4}\neq 0$, such that $\varepsilon_{4}=0.15$, and all
other $\tilde{\varepsilon}_{k}=0$, while the bottom panels are for a non-zero
$\tilde{\varepsilon}_{2}$ (with $\varepsilon_{2}=0.15$) and vanishing other
eccentricities. Figure 5 displays the time evolution of $v_{4}$ for systems
with $N_{\rm resc.}\approx 0.02$ (left) or 0.14 (right) rescatterings per
particle, while results for $N_{\rm resc.}\approx 0.35$ are shown in Fig. 9.

Overall, the results in Fig. 5 show that in the case of quadrangular flow
$v_{4}$, either from $\varepsilon_{2}$ or from $\varepsilon_{4}$, the $2\to 0$
model represents a very good approximation of the $2\to 2$ collision kernel
for low $N_{\rm resc.}$. In turn, the nice agreement with the analytical
results reinforce that statement and show that $v_{4}$ is proportional to
$\sigma$ in that regime. Indeed, the less good agreement of the “loss term”
results with the $2\to 0$ kernel for $N_{\rm resc.}=0.14$ can be attributed to
the limitation of the analytical calculations to linear order in the cross
section.

Although the results of Fig. 5 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_{\rm resc.}$ regime —, still there are important
differences. A first one, to which we shall come back in Sect. III.5, 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_{\rm resc.}\gtrsim 5$, with the $2\to 2$ collision kernel, since the $2\to
0$ scenario makes no sense in that case) the $v_{4}$ resulting from an initial
$\varepsilon_{4}>0$ is positive at late times, as found also e.g. in Ref. [31,
17],666Strictly speaking, in Ref. [31] a different initial profile was used,
namely Eq. (26). In turn, the results of Ref. [17] are for energy-weighted
quadrangular flow $v_{4}^{E}$, but we checked that it behaves like $v_{4}$ in
our setup. It seems that our small $N_{\rm resc.}$ regime is actually beyond
the low-opacity region studied in Ref. [17]. but contrary to the behavior of
the upper panels of Fig. 5. This means that the linear scaling with $N_{\rm
resc.}$ of the “final” $v_{4}$ observed in Fig. 5 breaks down at larger cross
sections. Note that a negative $v_{4}$ — more accurately, $v_{4}^{E}$ — for a
positive $\varepsilon_{4}$ in the few-rescatterings 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 times — namely a negative $v_{4}$ — in
case the system is initially deformed elliptically ($\varepsilon_{2}\neq 0$,
$\varepsilon_{4}=0$).

All in all, it seems that in the few-rescatterings regime $v_{4}$, either
resulting “linearly” from an initial quadrangularity $\varepsilon_{4}$ or
nonlinearly from an initial ellipticity $\varepsilon_{2}$, behaves like
elliptic flow $v_{2}$, i.e. it largely arises from the anisotropic escape of
particles. Interestingly, the contributions from $\varepsilon_{2}$ and
$\varepsilon_{4}$ to $v_{4}$ are of the same order of magnitude, and in the
small $N_{\rm 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.

###  III.4 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 $\sigma$. In turn, this
means that in the $2\to 0$ scenario $v_{5}$ arises at order $N_{\rm
resc.}^{2}$, while it is proportional to $N_{\rm resc.}$ in the $2\to 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.

Figure 6: Time dependence of hexagonal flow $v_{6}$ in systems with on average
$N_{\rm resc.}\approx 0.02$ per particle and with different initial
geometrical profiles, as described in the text.

Here we present in Fig. 6 results for $v_{6}$, which is at the limit of what
we can do numerically with reasonable control on the signal when $N_{\rm
resc.}\approx 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
an hexagonal deformation (linear response $v_{6}\propto\varepsilon_{6}$, dot-
dashed line), only an initial triangularity (quadratic response
$v_{6}\propto\varepsilon_{3}^{2}$, dashed lines), only an initial
$\varepsilon_{2}$ (cubic response $v_{6}\propto\varepsilon_{2}^{3}$, dotted
lines), or with both initial $\varepsilon_{2}$ and $\varepsilon_{4}$ (mixed
quadratic response $v_{6}\propto\varepsilon_{2}\varepsilon_{4}$, full line).
In every setup the only non-zero $\varepsilon_{n}$ are set to $0.15$. The
numerical results with an initial $\varepsilon_{2}$ resp. $\varepsilon_{3}$
are from the same simulations as in Sect. III.1 resp. III.2. We did not
attempt to perform simulations with an initial non-zero $\varepsilon_{6}$ nor
with both $\varepsilon_{2}$ and $\varepsilon_{4}$ (and aligned symmetry planes
$\Phi_{2}$ and $\Phi_{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 $\varepsilon_{3}$ agree rather well across the three
scenarios of this paper in the few-rescatterings regime. This agreement should
be contrasted with Sect. III.2, 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 $\varepsilon_{2}$, the results from
numerical simulations are extremely small but seem to be non-zero and
consistent in the $2\to 2$ and $2\to 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 $\varepsilon_{2}^{3}$ to $v_{6}$ (or to $v_{2}$) can only
arise at order $\sigma^{2}$, not at linear order in $\sigma$ as considered
here.

Eventually, the analytical results for initial geometries with either
$\varepsilon_{6}=0.15$ or $\varepsilon_{2}=\varepsilon_{4}=0.15$ are of the
same typical magnitude as those for $\varepsilon_{3}=0.15$, and as in the case
of $v_{4}$ the signal changes sign (here twice) over the evolution.

###  III.5 Local production rate of anisotropic flow

Figure 7: Angle-averaged local production rates $D_{2}$ (upper left), $D_{3}$
(upper right), and $D_{4}$ for an initial state with $\varepsilon_{2}=0$ and
$\varepsilon_{4}=0.15$ (bottom left) or with $\varepsilon_{2}=0.15$ and
$\varepsilon_{4}=0$ (bottom right).

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. (9), from which its expression can be read off. Figure 7
shows the results of our analytical approach for $D_{2}$, $D_{3}$, and
$D_{4}$, for the setups of Sects. III.1–III.3 with $N_{\rm resc.}=0.02$. Note
that we display $D_{n}(t,r)$ multiplied by $r$, so that the production rate
$\partial_{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 $\varepsilon_{n}$
exhibit similar qualitative features. The innermost region of the system —
extending up to $r\simeq R$ in the case of $D_{2}$, up to $r\simeq 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. 5, with a derivative that turns
negative around $t/R\simeq 1.5$, reflects the progressive dominance of the
region at intermediate $r$ in the lower left panel of Fig. 7. Similarly, the
(small) decrease of $v_{2}(t)$ for $t/R\gtrsim 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 $\partial_{t}v_{3}(t)=0$, while in Ref. [18] an “almost nearly perfect
cancellation” resulting in a very small negative $v_{3}^{E}$ 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)\propto t^{n+1}$ in the few collisions regime [32].

The lower right plot of Fig. 7, showing $D_{4}$ for the nonlinear response of
$v_{4}$ to an initial $\varepsilon_{2}$, is completely different, with a clear
negative contribution at early times and for $r\lesssim 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 $\varepsilon_{2}\neq 0$,
irrespective of its sign, i.e. the ellipse orientation, generates via the loss
term a negative $v_{4}\propto-\varepsilon_{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 $\varepsilon_{4}$, which is the seed
for the positive contributions to $v_{4}$ at later times, as seen in the
bottom panels of Fig. 5.

##  IV Summary

We have investigated anisotropic flow in three models: 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\to 0$ version of the code. Eventually, we carried out
analytical calculations within Boltzmann kinetic theory, including only the
loss term of the binary collision kernel and restricting ourselves to leading
order in the cross section.

The strength of the latter calculations is that they provide us directly with
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
yields $v_{n}=0$ for all odd coefficients, but not for even ones, which hints
at a fundamental difference between odd and even harmonics.

For even harmonics ($v_{2}$, $v_{4}$, $v_{6}$), the results of the three
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}\propto\varepsilon_{n}$ and the nonlinear response like e.g.
$v_{4}\propto\varepsilon_{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\to 2$ and $2\to
0$ numerical scenarios differ significantly, even in the very few
rescatterings regime. Indeed, the former scale roughly linearly with $N_{\rm
resc.}$, while the latter rather scale like $N_{\rm resc.}^{2}$. That finding
in the $2\to 0$ model is consistent with the fact that we find $v_{3}=0$ in
our analytical calculations at order $\sigma$ — while we found in a parallel
study that there is a non-zero $v_{3}$ at order $\sigma^{2}$ [32]. 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
probability, but that the fate of particles after they have undergone a
collision does matter.

We would thus conclude that the “escape mechanism” cannot explain the whole
anisotropic flow signal, but possibly only the even harmonics. 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,777A difference between even and odd harmonics of two-particle
azimuthal correlations was found in a study of proton–nucleus collisions
within a Color Glass Condensate based approach [39]. but may be relevant in
systems with small multiplicity. Our study within a toy transport model (two-
dimensional expansion, hard spheres) clearly needs to be replicated with more
realistic codes — which however have to be “truncated” to accommodate the
$2\to 0$ scenario.

###### Acknowledgements.

We thank Marc Borrell, Kai Gallmeister, Carsten Greiner, Sören Schlichting and
Clemens Werthmann for fruitful discussions. The authors acknowledge support by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through
the CRC-TR 211 ’Strong-interaction matter under extreme conditions’ - project
number 315477589 - TRR 211. Numerical simulations presented in this work were
performed at the Paderborn Center for Parallel Computing (PC2) and the
Bielefeld GPU Cluster, and we gratefully acknowledge their support.

##  Appendix A Odd flow harmonics in the “loss term” scenario

In this Appendix, we show that the odd flow harmonics $v_{3}$, $v_{5}$… are
identically zero when computed at leading order in the cross section with a
collision kernel including only the loss term of the Boltzmann equation, Eq.
(2).

Following Eq. (9), our analytical calculation of $v_{n}(t)$ involves the
integral over the transverse plane

${\cal I}\equiv\\!\int\\!f^{(0)\\!}({\bf x}-{\bf v}t,{\bf p})f^{(0)\\!}({\bf
x}-{\bf v}_{1}t,{\bf p}_{1})\,{\rm d}^{2}{\bf x},$ (20)

where the free-streaming distribution has been expressed in terms of the
initial condition via Eq. (6). A straightforward change of variable yields

${\cal I}=\\!\int\\!f^{(0)\\!}({\bf x}-\bm{\xi},{\bf p})f^{(0)\\!}({\bf
x}+\bm{\xi},{\bf p}_{1})\,{\rm d}^{2}{\bf x}$ (21)

with $\bm{\xi}\equiv\frac{1}{2}({\bf v}+{\bf v}_{1})t$.

Let us assume right away that the initial-state phase space density factorizes
into independent spatial and momentum distributions as in Eq. (10).888The
assumption only matters when going from Eq. (23) to Eq. (25). In that case the
momentum parts are irrelevant for the integral over ${\bf x}$ and the
calculation of ${\cal I}$ involves that of

${\cal I}^{\prime}=\\!\int\\!F({\bf x}-\bm{\xi})F({\bf x}+\bm{\xi})\,{\rm
d}^{2}{\bf x}.$ (22)

Since the integral runs over the whole transverse plane, we may equivalently
replace the integrand by its even part:

${\cal I}^{\prime}=\frac{1}{2}\\!\int\\!\big{[}F({\bf x}\\!-\\!\bm{\xi})F({\bf
x}\\!+\\!\bm{\xi})+F(-{\bf x}\\!-\\!\bm{\xi})F(-{\bf
x}\\!+\\!\bm{\xi})\big{]}{\rm d}^{2}{\bf x}.$ (23)

This integral can be further transformed by introducing the even and odd parts
of the spatial profile $F$:

$F({\bf x})=F_{+}({\bf x})+F_{-}({\bf x})\text{ with }F_{\pm}(-{\bf x})=\pm
F_{\pm}({\bf x}).$ (24)

The even “eccentricities” of the geometry and its isotropic component are
entirely controlled by $F_{+}$, while $F_{-}$ accounts for the odd
eccentricities. Replacing $F$ by $F_{+}+F_{-}$ in Eq. (23), the integrand
yields 8 terms: four of those cancel pairwise and there remains

${\cal I}^{\prime}=\\!\int\\!\big{[}F_{+}({\bf x}-\bm{\xi})F_{+}({\bf
x}+\bm{\xi})+F_{-}({\bf x}-\bm{\xi})F_{-}({\bf x}+\bm{\xi})\big{]}{\rm
d}^{2}{\bf x}.$ (25)

As is well established in model studies,999See also Ref. [32] for a more
formal proof within kinetic theory. in the absence of initial anisotropic flow
a given harmonic $v_{n}$ can only arise as linear response to a modulation of
the initial geometry in the same $n$-th harmonic — symbolically
$v_{n}\propto\varepsilon_{n}$ —, or as quadratic response to two geometrical
modulations that combine appropriately — symbolically
$v_{n}\propto\varepsilon_{k}\varepsilon_{n-k}$, or more generally (but this
case cannot be obtained in the analytical approach of the present paper)
$v_{n}\propto\varepsilon_{k_{1}}\cdots\varepsilon_{k_{m}}$ with
$k_{1}+\cdots+k_{m}=n$. Since all modulations of even (including 0) resp. odd
order are accounted for by $F_{+}$ resp. $F_{-}$, one sees that the products
$F_{+}F_{+}$ or $F_{-}F_{-}$ in the integrand of Eq. (25) can yield the
necessary contributions to $v_{n}$ of the kind $\varepsilon_{n}$ or
$\varepsilon_{k}\varepsilon_{n-k}$ for any even harmonic $n$, but not for odd
$n$. That is, the integral ${\cal I}^{\prime}$ does not depend on the momentum
azimuths $\varphi_{\bf p}$, $\varphi_{1}$ (on which $\bm{\xi}$ implicitly
depends) in such manner that after multiplying with the Møller velocity and
$\cos(n\varphi_{\bf p})$ and integrating over these azimuths, there could
result a non-zero $v_{n}$ when $n$ is odd.

To conclude, note that our proof does not explicitly use the dimensionality of
the system, nor does it make any assumption on the particle mass. However, it
assumes that the collision kernel does not include quantum effects.

##  Appendix B Results for $N_{\rm resc.}\approx 0.35$

Figure 8: Time dependence of elliptic flow $v_{2}$ in systems with initially
$\varepsilon_{2}=0.15$ and on average $N_{\rm resc.}\approx 0.35$
rescatterings per particle. The green curves are for systems with elastic
binary scatterings, the blue lines for the $2\to 0$ scenario and the red ones
show the analytical result (16).

Figure 9: Time-dependence of quadrangular flow $v_{4}$ in systems with $N_{\rm
resc.}\approx 0.35$ rescatterings per particle on average. Left: system with
an initial $\varepsilon_{4}=0.15$ and $\varepsilon_{2}=0$; right: system with
an initial state $\varepsilon_{2}=0.15$ and $\varepsilon_{4}=0$. The green
curves are for systems with elastic binary scatterings, the blue lines for the
$2\to 0$ scenario and the red ones show the analytical results (18) (left) or
Eq. (19) (right).

In this Appendix we provide for the sake of reference results for $v_{2}$
(Fig. 8) and $v_{4}$ (Fig. 9) for systems in which the mean number of
rescatterings per particle is about 0.35. For the $2\to 0$ scenario, this
means that approximately $70\%$ of the particles disappear over the system
evolution. Accordingly the assumption underlying the analytical calculations,
that the phase-space distribution $f(t,{\bf x},{\bf p})$ deviates negligibly
at all times from the free-streaming distribution $f_{\rm f.s.}(t,{\bf x},{\bf
p})$ with the same initial condition, is clearly non fulfilled.

The two plots displaying “linear” flow response, namely $v_{n}$ for an initial
non-zero $\varepsilon_{n}$ with $n=2$ (Fig. 8) or $n=4$ (Fig. 9 left), are
similar: The results from the simulations with the $2\to 2$ collision kernel
(green lines) and the $2\to 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\to 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.3,
the agreement for $v_{4}$ disappears at higher $N_{\rm resc.}$ values, since
the $2\to 2$ results become positive.

Going to the right panel of Fig. 9 showing the nonlinear response
$v_{4}\propto\varepsilon_{2}^{2}$, we just note that the agreement between the
three approaches is quite good, again without attempting to interpret it.

##  Appendix C Alternative distribution function

To check whether our main results are specific to our choice of initial
profile (11), we repeated the calculations in the few-rescatterings regime
$N_{\rm resc.}\approx 0.02$ for a slightly different geometry, namely

$\tilde{F}(r,\theta)=\frac{N\sqrt{1-\delta_{j}^{2}}}{2\pi
R^{2}}\exp\bigg{(}\\!\\!-\\!\frac{r^{2}[1+\delta_{j}\cos(j\theta)]}{2R^{2}}\bigg{)},$
(26)

which was used in the fluid-dynamical regime in Ref. [31]. A drawback of this
distribution is that a given $\delta_{j}$ contributes to several
eccentricities $\varepsilon_{n}$, namely for all $n$ that are multiples of
$j$. On the other hand, the density (26) is positive definite irrespective of
the parameter values.

As we could not perform all integrals with the distribution function (26)
analytically, we used the Vegas Monte Carlo (MC) integration method [40] to
evaluate the flow coefficients. In addition, we performed transport
simulations only with the $2\to 2$ collision kernel.

Figure 10: Time dependence of anisotropic flow harmonics for the initial
geometric profile (26) in systems with $N_{\rm resc.}\approx 0.02$. Left:
$v_{2}$ for $\delta_{2}=0.15$ (green: numerical simulations, red: semi-
analytical approach). Right: $v_{3}$ for $\delta_{3}=0.15$ (magenta), $v_{4}$
(orange: for $\delta_{2}=0.15$, purple: for $\delta_{4}=0.15$), and $v_{5}$
(cyan, for $\delta_{5}=0.15$) in numerical simulations (dashed lines) or
computed with Vegas (full lines; the color bands show the $3\sigma$ error of
the MC integration).

Figure 10 shows our results for the flow coefficients $v_{n}$ with
$n\in\\{2,3,4,5\\}$. The odd harmonics ($v_{3}$, $v_{5}$) from the semi-
analytical approach are zero within the error bars of the MC integration, as
expected from Appendix A, while the triangular flow $v_{3}$ is clearly non-
zero in the transport simulations. In contrast to this mismatch for the odd
harmonics, the $v_{2}$ and $v_{4}$ signals in a system with an initial
$\delta_{2}$ are in nice agreement in the two approaches. The agreement is
less good for the $v_{4}$ from an initial $\delta_{4}$, but this may be due to
the numerical noise in the simulations at early times, since at later times
the two curves run parallel to each other. In addition, the overall shapes of
$v_{2}$ and $v_{4}$ (either from a non-zero $\delta_{2}$ or a non-zero
$\delta_{4}$) are similar to those found in Sects. III.1 and III.3 with the
distribution (11).

## References

  * [1] U. Heinz and R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013). arXiv:1301.2826 [nucl-th]
  * [2] S. Voloshin and Y. Zhang, Flow study in relativistic nuclear collisions by Fourier expansion of azimuthal particle distributions, Z. Phys. C 70, 665 (1996). arXiv:hep-ph/9407282
  * [3] J. L. Nagle and W. A. Zajc, Small System Collectivity in Relativistic Hadronic and Nuclear Collisions, Ann. Rev. Nucl. Part. Sci. 68, 211 (2018). arXiv:1801.03477 [nucl-ex]
  * [4] S. Jeon and U. Heinz, Introduction to Hydrodynamics, Int. J. Mod. Phys. E 24, 1530010 (2015). arXiv:1503.03931 [hep-ph]
  * [5] P. Romatschke and U. Romatschke, Relativistic fluid dynamics in and out of equilibrium (Cambridge University Press, Cambridge, 2019). arXiv:1712.05815 [nucl-th]
  * [6] J. Y. Ollitrault, Anisotropy as a signature of transverse collective flow, Phys. Rev. D 46, 229 (1992).
  * [7] P. Huovinen and P. V. Ruuskanen, Hydrodynamic models for heavy ion collisions, Ann. Rev. Nucl. Part. Sci. 56, 163 (2006). arXiv:nucl-th/0605008
  * [8] D. A. Teaney, Viscous hydrodynamics and the Quark Gluon Plasma, in Quark Gluon Plasma 4, (editors R. C. Hwa and X. N. Wang, World Scientific, Singapore, 2010), p. 207. arXiv:0905.2433 [nucl-th]
  * [9] M. Luzum and H. Petersen, Initial state fluctuations and final state correlations in relativistic heavy-ion collisions, J. Phys. G 41, 063102 (2014). arXiv:1312.5503 [nucl-th]
  * [10] R. D. Weller and P. Romatschke, One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at $\sqrt{s}=5.02$ TeV, Phys. Lett. B 774, 351 (2017). arXiv:1701.07145 [nucl-th]
  * [11] W. Zhao, Y. Zhou, K. Murase and H. Song, Searching for small droplets of hydrodynamic fluid in proton–proton collisions at the LHC, Eur. Phys. J. C 80, 846 (2020). arXiv:2001.06742 [nucl-th]
  * [12] C. Gombeaud and J.-Y. Ollitrault, Elliptic flow in transport theory and hydrodynamics, Phys. Rev. C 77, 054904 (2008). arXiv:nucl-th/0702075
  * [13] P. Romatschke, Azimuthal anisotropies at high momentum from purely non-hydrodynamic transport, Eur. Phys. J. C 78, 636 (2018). arXiv:1802.06804 [nucl-th]
  * [14] N. Borghini, S. Feld and N. Kersting, Scaling behavior of anisotropic flow harmonics in the far from equilibrium regime, Eur. Phys. J. C 78, 832 (2018). arXiv:1804.05729 [nucl-th]
  * [15] A. Kurkela, U. A. Wiedemann and B. Wu, Opacity dependence of elliptic flow in kinetic theory, Eur. Phys. J. C 79, 759 (2019). arXiv:1805.04081 [hep-ph]
  * [16] A. Kurkela, U. A. Wiedemann and B. Wu, Flow in AA and pA as an interplay of fluid-like and non-fluid like excitations, Eur. Phys. J. C 79, 965 (2019). arXiv:1905.05139 [hep-ph]
  * [17] A. Kurkela, S. F. Taghavi, U. A. Wiedemann and B. Wu, Hydrodynamization in systems with detailed transverse profiles, Phys. Lett. B 811, 135901 (2020). arXiv:2007.06851 [hep-ph]
  * [18] A. Kurkela, A. Mazeliauskas and R. Törnkvist, Collective flow in single-hit QCD kinetic theory, JHEP 11, 216 (2021). arXiv:2104.08179 [hep-ph]
  * [19] V. E. Ambruş, S. Schlichting and C. Werthmann, Development of transverse flow at small and large opacities in conformal kinetic theory, Phys. Rev. D 105, 014031 (2022). arXiv:2109.03290 [hep-ph]
  * [20] N. Borghini, M. Borrell and H. Roch, Early time behavior of spatial and momentum anisotropies in kinetic theory across different Knudsen numbers, arXiv:2201.13294 [nucl-th]
  * [21] L. He, T. Edmonds, Z. W. Lin, F. Liu, D. Molnar and F. Wang, Anisotropic parton escape is the dominant source of azimuthal anisotropy in transport models, Phys. Lett. B 753, 506 (2016). arXiv:1502.05572 [nucl-th]
  * [22] H. Heiselberg and A. M. Levy, Elliptic flow and HBT in noncentral nuclear collisions, Phys. Rev. C 59, 2716 (1999). arXiv:nucl-th/9812034
  * [23] N. Borghini and C. Gombeaud, Anisotropic flow far from equilibrium, Eur. Phys. J. C 71, 1612 (2011). arXiv:1012.0899 [nucl-th]
  * [24] X. N. Wang and F. Yuan, Azimuthal asymmetry of $J/\psi$ suppression in noncentral heavy ion collisions, Phys. Lett. B 540, 62 (2002). arXiv:nucl-th/0202018
  * [25] N. Borghini and C. Gombeaud, Anisotropic collective flow of a Lorentz gas, J. Phys. G 38, 124172 (2011). arXiv:1106.6009 [nucl-th]
  * [26] N. Kersting, N. Borghini and S. Feld, Anisotropic flow in the few collisions regime: application to bottomonia, MDPI Proc. 10, 16 (2019). arXiv:1811.06195 [nucl-th]
  * [27] B. Alver and G. Roland, Collision geometry fluctuations and triangular flow in heavy-ion collisions, Phys. Rev. C 81, 054905 (2010). arXiv:1003.0194 [nucl-th] Erratum: Phys. Rev. C 82, 039903 (2010)
  * [28] D. Teaney and L. Yan, Triangularity and dipole asymmetry in heavy ion collisions, Phys. Rev. C 83, 064904 (2011). arXiv:1010.1876 [nucl-th]
  * [29] F. G. Gardim, F. Grassi, M. Luzum and J.-Y. Ollitrault, Mapping the hydrodynamic response to the initial geometry in heavy-ion collisions, Phys. Rev. C 85, 024908 (2012). arXiv:1111.6538 [nucl-th]
  * [30] H. Roch and N. Borghini, Fluctuations of anisotropic flow from the finite number of rescatterings in a two-dimensional massless transport model, Eur. Phys. J. C 81, 380 (2021). arXiv:2012.02138 [nucl-th]
  * [31] B. H. Alver, C. Gombeaud, M. Luzum and J. Y. Ollitrault, Triangular flow in hydrodynamics and transport theory, Phys. Rev. C 82, 034913 (2010). arXiv:1007.5469 [nucl-th]
  * [32] M. Borrell and N. Borghini, Early time behavior of spatial and momentum anisotropies in a kinetic approach to nuclear collisions, arXiv:2109.15218 [nucl-th]
  * [33] N. Borghini and J.-Y. Ollitrault, Momentum spectra, anisotropic flow, and ideal fluids, Phys. Lett. B 642, 227 (2006). arXiv:nucl-th/0506045
  * [34] D. Teaney and L. Yan, Non linearities in the harmonic spectrum of heavy ion collisions with ideal and viscous hydrodynamics, Phys. Rev. C 86, 044908 (2012). arXiv:1206.1905 [nucl-th]
  * [35] H. Niemi, G. S. Denicol, H. Holopainen and P. Huovinen, Event-by-event distributions of azimuthal asymmetries in ultrarelativistic heavy-ion collisions, Phys. Rev. C 87, 054901 (2013). arXiv:1212.1008 [nucl-th]
  * [36] L. V. Bravina et al., Is hexagonal flow $v_{6}$ just a superposition of elliptic $v_{2}$ and triangular $v_{3}$ flows?, Phys. Rev. C 89, 024909 (2014). arXiv:1311.0747 [hep-ph]
  * [37] J. Qian, U. W. Heinz and J. Liu, Mode-coupling effects in anisotropic flow in heavy-ion collisions, Phys. Rev. C 93, 064901 (2016). arXiv:1602.02813 [nucl-th]
  * [38] G. Giacalone, L. Yan and J.-Y. Ollitrault, Nonlinear coupling of flow harmonics: Hexagonal flow and beyond, Phys. Rev. C 97, 054905 (2018). arXiv:1803.00253 [nucl-th]
  * [39] M. Mace, V. V. Skokov, P. Tribedy and R. Venugopalan, Systematics of azimuthal anisotropy harmonics in proton–nucleus collisions at the LHC from the Color Glass Condensate, Phys. Lett. B 788, 161 (2019). arXiv:1807.00825 [hep-ph] Erratum: Phys. Lett. B 799, 135006 (2019)
  * [40] P. Lepage, gplepage/vegas: vegas version 3.5.3, Zenodo (2020), doi:10.5281/zenodo.4268682