Collective flow in single-hit QCD kinetic theory

Motivated by recent interest in collectivity in small systems, we calculate the harmonic flow response to initial geometry deformations within weakly coupled QCD kinetic theory using the first correction to the free-streaming background. We derive a parametric scaling formula that relates harmonic flow in systems of different sizes and different generic initial gluon distributions. We comment on similarities and differences between the full QCD effective kinetic theory and the toy models used previously. Finally we calculate the centrality dependence of the integrated elliptic flow $v_2$ in oxygen-oxygen, proton-lead and proton-proton collision systems.


Introduction
The physical picture implemented in event generators used in phenomenological modelling of soft particle production in proton-proton collisions [1,2] is maximally different from that implemented in the descriptions of nucleus-nucleus collisions [3,4]. In the former, partons fragment but free-stream to the detector after being created in an initial hard scattering without undergoing secondary rescatterings with the other partons. In the latter, reinteractions are so strong that a nearly ideally hydrodynamized fluid is created. The recent experimental characterization of the smooth onset of collectivity, i.e., the onset of qualitative features characteristic to the fluid-dynamic picture and absent in the free-streaming picture, challenges this dichotomy of two disconnected pictures [5]. Prime examples of the signs of collectivity include the strangeness enhancement observed in proton-proton (pp), proton-nucleus (pA), and in nucleus-nucleus (AA) collisions [6] as well as the formation of multi-particle collective flow, measured by the azimuthal harmonics v n and their cumulants [7][8][9]. These observations offer an inroad to study how the ideal fluid is built up from fundamental interactions of elementary particles and how hydrodynamization takes place as a function of the system size [5].
There have been multiple attempts to describe collectivity in small systems extending the models of pp and AA collisions [10,11]. However, in order to fully exploit the experimental progress, the experimental data needs to be confronted with models that encompass the both extremes and are able to dynamically describe the process of hydrodynamization. One such a model is the QCD effective kinetic theory (EKT) [12] that gives a leading-order accurate description of the matter created in the asymptotic limit of infinite collision energy √ s → ∞. It builds on a picture arising from perturbative QCD in which partons undergoing elastic scattering and radiate medium-induced collinear radiation. This effective theory has been extensively employed to study how non-abelian gauge theories thermalize and how the hydrodynamization takes place in heavy-ion collisions in the weak coupling limit [13][14][15][16][17][18][19][20][21] (see also [22,23]). It has also been used to study how hard partons are hydrodynamized in the QCD medium leading to jet quenching [24,25] (see also [26][27][28][29]). Here, we will take the first steps to employ the EKT to study how hydrodynamization takes place as a function of system size and specifically how signals of collectivity arises in the smallest systems.
In the case of sufficiently small systems, the system stays far from equilibrium throughout the evolution and the deviation from collisionless expansion can be treated as a perturbation to the free-streaming background. This single-hit approximation has been studied for multiple simplified models of QCD medium capturing different aspects of the full effective kinetic theory, including the isotropization time approximation (ITA) [30][31][32][33] and other generic models of elastic scattering [34][35][36][36][37][38][39]. These works have demonstrated that the first scatterings have the largest effect in formation of the flow harmonics and that already the single-hit approximation leads to sizeable elliptic flow v 2 . This qualitative observation suggests that collectivity in a form of azimuthal flow is a signal of final state interactions rather than of full hydrodynamization. In ref. [33], the systematics of how different flow harmonics arise in kinetic theories have been characterized and contrasted with how they arise in relativistic fluid dynamics.
The QCD effective kinetic theory has more structure than the simple ITA kinetic theory. In ITA, the response of the flow coefficients to initial geometry depends only on a single parameter, opacityγ, that describes the system size in units of the mean free path [30]. The situation in EKT is somewhat more complicated because the collision kernel depends non-trivially on the interplay of classical and Bose-enhanced parts of the scattering terms as well as on the in-medium screening that regulates the soft divergence of the elastic scattering. Hence, in addition to depending on the system size, the flow is sensitive to the occupancy determining the strength of Bose enhancement as well as the coupling constant, determining the amount of Debye screening. Here, we will fully characterize the energy weighted linear flow response coefficient v 2 /ε 2 in terms of these variables.
We note that while EKT gives a leading order description of time-evolution of large and isotropic systems, in the current work we are pushing it beyond its strict applicability. Hence, our results are not fully leading-order accurate but represent rather a model calculation that is based on our current best understanding of perturbation theory. The manuscript is organized as follows. In Section 2 we summarize the model setup and discuss the computation of the linear response function v n /ε n using the co-moving coordinate system. In Section 3 we discuss the conformal scaling of EKT results and make comparisons to ITA. Finally, we present an exploratory study of elliptic flow signal in oxygen-oxygen (OO), proton-lead (pPb) and proton-proton (pp) collisions using single hit EKT.

Single hit approximation
We consider the leading order QCD effective kinetic theory formulation by Arnold, Moore and Yaffe [12]. In a boost invariant system the color and polarization averaged distribution f is governed by the Boltzmann equation [40] ∂ ∂τ is a sum of a leading-order elastic 2 ↔ 2 scattering kernel containing the physics of in-medium screening and an effective mediuminduced collinear 1 ↔ 2 radiation kernel. For the explicit discussion of the implementation we refer to previous papers [14,41]. In the following we will restrict ourselves only to a purely gluonic plasma. We are interested in small, dilute systems, where partons experience only a few scatterings. Therefore it is useful to first consider the solution to the collisionless Boltzmann equation f (0) . The free-streaming solution can be written entirely in terms of the initial distribution given at time τ 0 using a co-moving coordinate system [31] where the longitudinal momentum is simply the rescaledp z = p z τ τ 0 , while the transverse co-moving coordinate is In co-moving coordinates, defining the co-moving distribution functionf (τ,x ⊥ ; p ⊥ ,p z ) = f (τ, x ⊥ ; p ⊥ , p z ), the Boltzmann equation reduces to an ordinary differential equation in time τ , This form is convenient for linearizing the solution in the number of scatteringsf =f (0) + f (1) + . . ., where the background and the single-scattering distributions evolve according to This is the so-called single hit approximation and the first correction to the distribution can be obtained directly by integratioñ where the initial condition for the perturbation isf (1) (τ 0 ) = 0 by definition.

Energy flow
We will calculate the elliptic flow generated in the single-hit approximation in QCD kinetic theory. The elliptic flow v 2 is customarily defined by the azimuthal deformation of the particle number. The translation of partonic degrees of freedom to observed hadronic ones is sensitive to the details of the hadronization procedure as there is no conservation of particle number between the partonic and hadronic degrees of freedom of QCD. We instead consider the soft and collinear safe transverse energy flow which is insensitive to the hadronization effects. The transverse energy flow per unit rapidity reads as a function of the momentum-space azimuthal angle φ p , where p ⊥ = p ⊥ (cos(φ p ), sin(φ p )).
The energy weighted flow coefficients v n here are defined as Fourier components of the transverse energy flow where ψ n defines the azimuthal orientation of nth harmonic flow. Inserting the single hit expansion f ≈ f (0) + f (1) in Eq. (2.6) and noting that the symmetric free-streaming background has vanishing harmonics (assuming no initial momentumspace anisotropy), we get where we changed to co-movingx ⊥ ,p z coordinates and gained a τ 0 /τ factor as a Jacobian. Now we can substitute the single-hit approximation, Eq. (2.5), and revert to labcoordinates.
In leading order QCD the collision kernel consists of two terms: C = C 2↔2 + C 1↔2 . However the collinear 1 ↔ 2 splitting does not change the particle orientation and only shifts particles in energy. Therefore in the single hit approximation only 2 ↔ 2 scatterings will contribute to the soft and collinearly safe transverse energy weighted flow of Eq. (2.10) and momentum integrated v n is given by (2.11)

Linearized perturbation
We further simplify the problem by considering small azimuthal perturbations δf of an otherwise symmetric backgroundf , ). Then the collision kernel can be also linearized At initial time the background distributionf is azimuthally symmetric both in momentum and coordinate space, therefore scatterings of the background term C 2↔2 f will not contribute to anisotropy generation. Then, to linear order in perturbations, the momentum integrated v n is v n = − τ τ 0 dτ τ rdrD n (τ , r) (2.14) where we defined the partial integral D n as The linearized collision kernel δC 2↔2 f , δf can be obtained by straightforward linearization of the collision matrix |M| 2 (via m 2 g ) and the phase-space terms. For simple single harmonic perturbations of the background discussed below, the φ x integral in Eq. (2.15) can be performed explicitly. Then D n (τ, r) can be calculated by Monte-Carlo sampling of the multidimensional phase-space integral on a discrete τ − r grid. The integrated v n is then found by numerically integrating D n (τ, r). We refer to Appendix A for further details.

Elastic scattering kernel and screening mass
The elastic collision kernel is given by the multi-dimensional phase-space integral [12] C where the leading order scattering-amplitude-square is given in Mandelstam variables s, t, u and the 't Hooft coupling constant λ = N c g 2 as Small angle scattering then exchange momentum q = |p − p | → 0 and the t-channel (and similarly u-channel) suffers from Coloumb divergence. Medium screening effects have to be taken into account. This can be achieved by introducing a regularized t reg in the denominator of Eq. (2.17) (identically for u reg ), . (2.18) Here m 2 g is the gluon screening or effective mass The constant e 5/3 /4 in Eq. (2.18) is chosen to reproduce gluon drag and diffusion in full HTL treatment near equilibrium [13]. We note that the effective mass is Lorentz invariant, which follows from the fact that d 3 p/|p| is a Lorentz invariant integration measure and the phase space density is a Lorentz scalar [42]. Therefore we can perform the integral in Eq. (2.19) in the lab frame.

Background
For initial conditions of the spatially azimuthally symmetric background at an initial time τ 0 we consider the following ansatz designed to capture certain properties of the gluon distribution in a gluon-saturation, or Color-Glass-Condensate (CGC), -based calculation [14,43], where A controls the magnitude of the occupation. In CGC calculations in the weakcoupling limit A ∝ 1/λ 1. The parameter ξ controls the longitudinal momentum asymmetry, which is assumed to be large in CGC-type initial conditions. Here, Q(x ⊥ ) is the characteristic energy scale, which can be related to the average transverse momentum (2.21) For the transverse Q profile we take a simple Gaussian of width √ 2R 0 , Defining a short-hand notationÂ = Aτ 0 ξR 0 , the initial (lab-frame) energy density per rapidity (for ξ 1) is while the transverse energy integral (conserved by free-streaming) is Here ν g is the number of gluon degrees of freedom (ν g = 16 for N c = 3).
The elastic scattering rate depends on the gluon screening mass, Eq. (2.19). For an anisotropic distribution ξ 1 and at the initial time τ 0 , it can be written as To study the sensitivity of our results to the precise form of the distribution function, we will also consider a deformed Bose-Einstein distribution where 8ζ(5)/ζ(3) ≈ 2.6 preserves the equality between Q and the transverse momentum at the initial time.

Perturbation
For the linearized azimuthal perturbation we take where 1 is a number characterising the size of the perturbation (and which can be scaled out from calculations in linearized equations). Note that such a choice of perturbation in Eq. (2.27) implies that, by symmetry, ψ n = 0 in Eqs. (2.7) and (2.15) for the orientation of the flow harmonics.
We will mainly present results for the elliptic deformation n = 2. For such a perturbation we can define the eccentricity with respect to the background energy density which has an opposite sign to . As geometry anisotropy is transformed to momentum anisotropy with a sign change, such definition of geometry eccentricty keeps the response coefficient v 2 /ε 2 positive. Similarly we can define n = 3 eccentricity as (2.29)

Conformal scaling of the solutions
The EKT of Eq. (2.1) possesses conformal symmetry which connects the flow from simulations performed with different parameters. The different parameters characterizing the initial condition are the spatial extent (R 0 ), mean transverse momentum (Q 0 ), and initial occupancy (A). In addition, the initial condition depends on the initial anisotropy (ξ) at the initialization time (τ 0 ). We note that the kinetic theory itself depends on the coupling λ, both directly setting the overall rate of the collisions as well as determining the ratio m g /Q 0 . Lastly, to linear order in perturbations the flow response is linearly proportional to initial geometric deformations characterized by eccentricities ε n . We first discuss how simulations with different initial conditions are connected to each other through conformal scaling followed by a discussion about how varying the coupling of the kinetic theory affects the flow variables. The conformal scaling and the lack of intrinsic scales allows to trivially relate all simulations with different spatial extents and initial transverse momentum magnitude. This can be seen by scaling all coordinate and momentum space variables in Eqs. (2.14) and (2.15) by appropriate powers of R 0 and Q 0 respectively, i.e.,τ = τ /R 0 andr = r/R 0 , so that where the rescaledD n iŝ Here,Ê ⊥ = E ⊥ /R 3 0 Q 4 0 is the dimensionless transverse energy. Recall that the collision kernel has units of momentum and δĈ = δC/Q 0 has absorbed one power of Q 0 . After rescaling the dimensionful quantities,v n = v n /(R 0 Q 0 ) depends only on dimensionless parameters, A, ξ andτ 0 = τ 0 /R 0 .
If we start with the initial particle distribution isotropic in the transverse momentum space, but anisotropic in coordinate space, most of the final v n will be generated at times τ ∼ R 0 when particles from different density regions will rescatter. In particular, we do not expect flow generation at the initial time τ 0 R 0 and we would therefore like to take thê τ 0 → 0 limit. At early timesτ 1 when the transverse expansion may be neglected, the main effect of the free-streaming expansion is the linear increase of longitudinal anisotropy with time. Therefore a well-behaved limit τ 0 → 0 can be taken keepingτ 0 /ξ fixed.
Furthermore, we note that if the momentum space distribution is anisotropic enough, ξ 1, it essentially behaves as a delta function in p z , f ∝ δ(p z ). In this case the dependence in the overall occupancy of the initial distribution, A, can be absorbed to either initial timê τ 0 or equivalently to initial anisotropy ξ. That is, if there are two initial conditions with a fixed A/ξ but differing A and ξ, the one with the larger A will have a larger occupancy but less phase space occupied so that both distributions have the same number of particles and other integral moments. Therefore, keeping the parameters of the theory fixed (i.e., λ) but varying the initial conditions, the rescaled flow variablev n depends only on the scaling variableÂ assuming that τ 0 R 0 and ξ 1. Here we discuss the relation of the scaling variableÂ to initial gluon multiplicity for CGC-type initial conditions. In the saturation framework, the initial gluon multiplicity per rapidity and unit transverse area is given by [44,45] where c ≈ O(1) is gluon liberation coefficient [43]. On the other hand, the total particle number integral in our parametrization scales as ThereforeÂ has the following dependence on initial gluon multiplicity for CGC-type initial conditionsÂ We note that perturbatively 1/λ ∝ log (Q 0 /Λ), where Λ is the non-perturbative QCD scale.
For a system composed of nucleons R 0 > 1/Λ and in the perturbative Q 0 → ∞ limit we therefore get thatÂ 1. Now we consider the scaling of the numerator in Eq. (3.2). The loss and gain terms in the collision integral consist of a classical piece, which is quadratic in phase-space distributions and a cubic Bose-enhanced piece. These two terms scale as A 2 and A 3 , respectively. As the dependence on the initial conditions arises only through scaling variables, the dependence on A implies that the two terms are proportional toÂ 2 andÂ 3 . Cancelling one power ofÂ in the denominator, we can express the flow v n as a sum of classical and Bose-enhanced terms: The above Eq. (3.8) describes how the flow variables change when the initial conditions are varied but the coupling of the kinetic theory itself is kept fixed. Now we will finally discuss how different simulations with different values of λ are connected to each other. First we note that the matrix element is explicitly proportional to λ 2 so thatṽ 2 ≡ λ 2v 2 ∼ λ 2 . In addition to this explicit dependence, the in-medium screening regulating the matrix element depends on the coupling: a more strongly coupled medium leads to more screening and hence less small-angle scattering. This dependence on m 2 g /Q 2 0 is genuinely non-trivial   and simulations with different m 2 g /Q 2 0 cannot be related to each other through a simple scaling. Therefore bothv cl. n andv b.e.
n are non-trivial functions of this ratio. As the dependence on the initial conditions must come in the specific combination set by the scaling variable of Eq. (3.3), the dependence on m 2 g /Q 2 0 can only enter througĥ where we used the fact that m 2 g (τ 0 , r = 0) ∼ λAQ 2 0 /ξ, see Eq. (2.25). We finally note that the two scaling variablesÂ andm 2 g are related to each other for a given coupling λ as the amount of screening is related to the occupancyÂ. For the distribution given by λÂ, (3.10) as is evident from Eq. (2.25).
Our final scaling formula for the linear flow response in single-hit EKT, depends on two unknown functionsv cl. n (m 2 g ) andv b.e. n (m 2 g ). In addition, we identify the overall factorR = λ 2 R 0 Q 0Â as a parametric estimate of the system size in units of the (transport, large-angle) mean free path at the origin of the system, r = 0, at the time τ = R 0 when the flow is built up, where ξR 0 /τ 0 is the anisoptropy of the distribution at time τ = R 0 . Alternatively, we can understand Eq. (3.11) in terms of initial gluon multiplicity. For the CGC-type initial conditions (see Eq. (3.6)) and constant λ, we have thatm 2 g ∝Â ∝ (dN g /dη) −1/2 , whileR is approximately independent of multiplicity.
In Fig. 1 we display the classicalv cl. n (m 2 g ) and the bose-enhancedv b.e. n (m 2 g ) scaling functions as a function ofm 2 g . We observe that in both left and right panels different simulations collapse to a single universal curve within the statistical uncertainties. This demonstrates that the scaling form of Eq. (3.11) derived in τ 0 R 0 and ξ 1 is born out of realistic EKT simulations. This is a highly non-trivial validation of the scaling formula. In order to obtain a pocket formula for the EKT response, we perform a linear fit in the loglog plot, which reasonably well describes the scaling functions. We note that for parameter choices far from this scaling regime (see Table 3), the data points do not fall on the same universal curve. Such points outside of range of the scaling regime (labelled as OR for "Outside of Range") are shown in gray in Fig. 1, but are not included in the fit. We note that although the scaled classical part in Fig. 1 is numerically smaller than the Bose-enhanced term, the latter is normally suppressed byÂ 1 in Eq. (3.11). For example, for λ = 10, one needsm 2 g > 0.6 to getÂ > 1.

Comparison to Isotropization Time Approximation
In the Relaxation Time Approximation (RTA) the collision kernel is taken to be [46] − C[f ] = −γe where γ is the single parameter controlling the rate of relaxation, e rest is the (rest-frame) energy density in the frame moving with velocity u µ , which solves the Landau matching condition: T µν u ν = −e rest u µ , andĒ p = −p µ u µ is the particle energy in that frame. The term f iso (Ē p ; ε) is an isotropic thermal distribution towards which the system is relaxing, e.g., Boltzmann or Bose-Einstein distributions. The evolution of the energy flow in RTA does not depend on the p-dependence of f iso , but only on the energy density, rest-frame and the assumption that f iso is isotropic. Therefore it is sufficient to study the momentum integrated Boltzmann equation and only keep track of the isotropization in the Isotropization Time Approximation (ITA) [30]. In this section we bring our EKT results in contact with previous studies of 2D kinetic theory using ITA kinetic theory [20,[30][31][32][33]. In ref. [30] it was noted that the ITA flow response is uniquely determined by the opacity parameterγ ≡ γ(e rest τ )   units of the mean-free-path, c.f.R in Eq. (3.12). A fair comparison between EKT and RTA should then be done for the same mean-free-path length. However, it is not trivial to relate the normalization of l mfp in the two kinetic theories and for simplicity we will compare ITA and EKT at fixedγ. To make a connection with EKT, we replace γ with the specific shear viscosity η/s using the equilibrium relation in ITA kinetic theory: γ = T In the linear regime, the ITA response for the integrated harmonic flow is given by the coefficients [30] v ITA We have implemented the RTA collision kernel, Eq. (3.13), and verified that we reproduce the results above obtained by the momentum-integrated Boltzmann equation in ITA.
The collision kernel of RTA in Eq. (3.13) is much simplier than the EKT collision kernel, Eq. (2.16). In particular, the linear flow response dependence on the occupation is onlyÂ 1/4 , see Eq. (3.14), in contrast to Eq. (3.11), which depends non-trivially on m 2 g . In Fig. 2 we show the elliptic flow response v 2 /(ε 2γ ) in EKT as a function ofm g . The points correspond to EKT simulations for λ = 10. The black lines are the results obtained by summing the power-law fits in Fig. 1 for different values of λ = 5, 10, 20 with η/s ≈ 2, 0.6, 0.2. For reference, we display the ITA value by a horizontal dot-dashed line. The ratio with the opacityγ partially, but not completely, cancels the strong dependence A simple cartoon illustrating the pertinent features of elliptic and triangular flow generation in the single-hit approximation. The initial anisotropy can be crudely modelled by 2 or 3 hot spots (blue stars), and the free streaming particles can be followed; the dashed circles show the location of the particles at some later time. For v 2 , a majority of the flow is generated around the origin (red point), where the local isotropisation of the right and left movers leads to a global excess of up and down movers. As the momentum space anisotropy at the origin is aligned with the initial spatial anisotropy, the resulting momentum space anisotropy is anti-aligned with the spatial anisotropy. For v 3 , in addition to the origin, particles scatter also at r > R 0 at locations marked by the red dots. The local anisotropy in the origin is anti-aligned with the initial spatial anisotropy leading this time to an aligned momentum space anisotropy. At the red dots, the rest frame is moving away from the collision zone. The collisions at red dots lead to an excess of particles moving along with the rest frame leading to an anti-aligned contribution. For v 3 the two regions (green dot, red dots) approximately cancel each other, leading to a small total contribution. on the coupling constant λ in Eq. (3.11). The EKT result for v 2 /(ε 2γ ) grows approximately linearly withm g and form g ≈ 0.15 it is roughly equal to the corresponding ITA values. At small values of λ or smallÂ values the EKT can be less efficient in generating elliptic flow than ITA at the sameγ value. In the next section we will see that for initial conditions found in small collision systems, the EKT response is similar to that in ITA.
We close this section with a comment on the triangular flow. All of the discussion for linear elliptic flow response also generalizes to arbitrary v n harmonic and, in particular, the triangular flow n = 3. However, numerically, it is more difficult to study the scaling properties of v 3 response due to large numerical cancellations in the integral of D n (τ, r) in Eq. (2.14). A heuristic argument for the difference between the elliptic and triangular flow generation is given in the caption of Fig. 3 for a toy example with n = 2 or n = 3 point sources.
In Fig. 4 we display the results for the differential v n /(ε nγ ) distribution, i.e. −rτD n (τ ,r)/(ε nγ ), for n = 2, 3 harmonics in ITA (top row) and EKT (bottom row). The integral of these distributions are equal to v 2 /(ε 2γ ) and v 3 /(ε 3γ ) correspondingly. Indeed, we observe that  for v 2 the distribution peaks at r ≈ τ ≈ R 0 with positive contributions for all r and τ values both for ITA and EKT. In contrast, the negative triangular flow response v 3 /ε 3 is generated at small radii r R 0 and positive only for r R 0 . Therefore there are significant cancellations between the two regions and the net response is small. In ITA the positive component is dominant and we obtained a net positive v 3 /(ε 3γ ). However, in EKT, the contributions from large radii are smaller and nearly perfectly cancel the negative component at small radii. The end result is that v 3 /(ε 3γ ) has a very small negative value. Finally, for completeness, in Fig. 5 we show triangular flow results for a particular set of CGC-type initial conditions. We see that both the classical and Bose-enhanced parts are an order of magnitude smaller than the corresponding terms for the elliptic flow.

Energy weighted elliptic flow in small systems from single-hit EKT
In this section we apply the scaling laws for EKT flow response extracted in Section 3.1 to realistic situations that take place in ultra-relativistic collisions of light nuclei, protonnucleus and proton-proton collisions. In order to do so, we will use successful initial state, equilibration and hydrodynamic models to determine realistic initial conditions in small collision systems. We will then ask, what would be the expected v 2 signal if the system with the same initial conditions was to evolve in the single-hit EKT approximation. We emphasize that in this work we do not attempt to provide a complete description of signals of collectivity observed in small systems, as it clearly requires a detailed study of multiple observables. Rather, this is a proof-of-principle study of how efficient single-hit EKT is in generating elliptic flow signals.   The dynamical response of the EKT for a set of scaling variables is fully described by the scaling formula Eq. (3.11). The combined power-fit functions results in the following pocket formula for the elliptic flow 16) and the flow harmonic is given by the initial eccentricity ε 2 , dimensionless system sizeR, and dimensionlessÂ andm 2 g . These scaling parameters can be related to the physical dimensionful parameters (R 0 , (eτ ) 0 , and Q 0 ) via equations in Section 2.5.
We determine R 0 and (eτ ) 0 for different collision systems from the tabulated values of eccentricity ε 2 , RMS radius r 2 and entropy density dS/dy/(π r 2 ) (in arbitrary units) from ref. [47], that were generated using the TRENTo initial state model [48,49]. We summarize the procedure of obtaining initial conditions and tabulate the parameter values in Appendix C.
The knowledge of R 0 and (eτ ) 0 values is not, however, sufficient to uniquely determine Q 0 andÂ. This is because (eτ ) 0 ∝ÂQ 4 0 R 0 depends only on the combinationÂQ 4 0 and to break this degeneracy we need to also provide an estimate for the typical momentum scale Q 0 . In the saturation framework [44,45], the initial gluon multiplicity per unity rapidity and area is ∝ Q 2 0 (see Eq. (3.4)). The mean transverse momentum of such gluons is of the order of the saturation scale Q 0 . Therefore the initial gluon energy density in CGC-type initial conditions is We fix the proportionality coefficient in Eq. (

3.19)
These relations yield decreasing Q 0 , but increasingÂ in smaller collision systems. In Fig. 6 we illustrate the single-hit EKT response for elliptic flow in √ s NN = 5.02 TeV OO, pPb and pp collisions. These are small collision systems for which a single-scattering approximation might be more appropriate than an infinite rescattering limit found in an ideal hydrodynamic description. In the left panel we show the centrality dependence of v 2 in different systems for λ = 10. Q 0 decreases from ≈ 1.5 GeV in central to ≈ 0.3 GeV in peripheral collisions for all three systems (see Tables 5, 6 and 7). CorrespondinglyÂ increases from ≈ 0.02 to ≈ 0.2. The single-hit EKT response per eccentricity and system size, i.e. v 2 /(ε 2 R 0 Q 0 ), is a function ofm 2 g ∝Â, with stronger response for largerÂ (see Fig. 2). Even though in more peripheral bins for OO the system size and Q 0 reduces, the increase inÂ and ε 2 results in a weak centrality dependence of net v 2 . For comparison we also show the ITA response, Eq. (3.15), for the same η/s = 0.6 value. In the ITA case v 2 /(ε 2γ ) is constant and therefore the v 2 becomes small in peripheral bins. We find that for λ = 10 and Q 0 = 3 GeV the EKT and ITA response are similar in magnitude in central OO collisions, but deviate in peripheral bins.
In the right panel of Fig. 6 we study the sensitivity of the EKT response to the coupling constant λ and scale Q 0 in OO collisions. Green, blue and red lines correspond to λ = 20, 10, 5 respectively and we see a strong dependence on the coupling constant. For λ = 10, we also show the blue band obtained by varying Q 0 by 33% (i.e. varying Q 0 between 2 GeV and 4 GeV in a central PbPb collision). It is clear that the λ and Q 0 values in elliptic flow response are degenerate. In addition we show the centrality dependence of the initial geometry eccentricity scaled by the ideal (and conformal) hydrodynamic response [33]. This illustrates the diametrically opposite limit of infinite rescatterings. As flow then follows the eccentricity, we observe that it is larger in more peripheral bins and is generally larger than the single hit EKT response for given values of λ and Q 0 .
Experimentally measured charged particle number elliptic flow in pPb and pp collisions is only weakly dependent on centrality and of typical size v 2 {2} ≈ 0.06 [50]. We see in Fig. 6 that for our choice of input parameters, the EKT single-hit approximation reproduces the same order of magnitude of the elliptic energy flow. Clearly, one should not expect that our simplified model would accurately describe the experimental data. It only demonstrates that approaches based on QCD effective kinetic theory can efficiently generate sizable harmonic flow even in the limit of few rescatterings.

Conclusions and Outlook
In this work we presented the first study of system-size dependence of harmonic flow response in QCD effective kinetic theory. We used the single-hit approximation to calculate the linear response coefficient for energy weighted elliptic flow on top of a free-streaming background. Despite the simplifying assumptions, our study addresses a number of new dynamical features of the system that were not accessible in previous toy models. Energy flow response in QCD EKT is generated via the elastic 2 ↔ 2 scatterings, which can be Bose-enhanced if the phase-space occupation density is large enough. In addition, the elastic scattering matrix element is regulated by the in-medium screening mass. This leads to a non-trivial scattering rate and flow response dependences on the initial conditions. We find the scaling laws that relate flow response with different initial conditions to each other and we provide a simple pocket formula parametrization of the elliptic flow response in EKT.
We apply single-hit EKT response to estimate the centrality dependence of energy weighted elliptic flow in √ s NN = 5.02 TeV OO, pPb and pp collisions. Although there are significant systematic uncertainties and simplifications involved, the resulting energy weighted elliptic flow for a realistic choice of parameters was found to be in order of magnitude agreement with the experimentally measured charged particle number elliptic flow in pPb and pp collisions. For the initial conditions studied in this paper both the EKT and ITA responses are rather similar for the energy weighted elliptic flow, but we found a much smaller energy weighted triangular flow in the EKT simulations than in those using ITA. This indicates that further study of kinetic theories with more complicated collision kernels like EKT can lead to specific collective flow signatures that may allow to distinguish between different microscopic interaction mechanisms. One of the clear advantages of EKT over ITA, is that the momentum dependence of flow harmonics can be studied. However, to correctly describe the p T -resolved v n one will need to take into account the p T profile of initial conditions, the collinear processes and hadronization. Finally, going beyond the linearized single-hit approximation employed in this work will be important for connecting the small and large systems and studying the hydrodynamization as a function of the system size.
[54] G. Giacalone, A. Mazeliauskas and S. Schlichting, Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions, Phys. A Linearizing D n In this section we derive the explicit expression for the time-and radius-resolved flow distribution D n = − dvn τ dτ rdr given in Eq. (2.15). The linearized elastic collision kernel δC 2↔2 contains two types of terms. The first of them is due to the linearization of the phase-space distributions in the loss and gain factors in Eq. (2.16). Using the free-streaming solution, we can write for time τ ≥ τ 0 Here the |x ⊥ | and φx is the radius and angle of the co-moving coordinate, Eq. (2.2). The second term in δC 2↔2 arises due to linear variation of the gluon screening mass in the regulated t reg and u reg , see Eq. (2.18), in the matrix element Eq. (2.17). Splitting the screening mass in the background and perturbation m 2 g =m 2 g + δm 2 g and using Eq. (2.19) we can write that In the second line we expanded the cos(nφx) = cos(nφ x ) cos(nφ ∆x ) − sin(nφ x ) sin(φ ∆x ), where φ ∆x = φx − φ x . We dropped the sin terms, because |x ⊥ | is an even function in the relative momentum angle φ p − φ x , while φ ∆x is odd. Finally in the last line we defined δm 2 g and explicitly factored out the angular dependence on φ x . Then we can linearize the scattering matrix as With the definitions given above the first step of the linearization of D n is now straight- We now turn our attention to the angular dependence of the equation above. Both terms of the distribution,f and δf , have an angular dependence through the magnitude of the co-moving coordinate |x ⊥ |. By squaring Eq. (2.2) we can write Hence, |x ⊥ | depends only on the relative angle φ p −φ x . The only other angular dependence appearing in Eq. (A.4) are the explicit terms cos(nφx (k) ) cos(nφ p ). Using again φ ∆x = φx − φ x , which is a function of the relative angle φ k − φ x , and writing φ p = φ p − φ x + φ x we factor cos(nφx) cos(nφ p ) in terms depending on the relative momentum angles and φ x as follows cos (nφx) cos(nφ p ) = = cos (nφ ∆x ) cos(n(φ p − φ x )) cos 2 (nφ x ) Hence, every term in the integral in Eq. (A.4) depends only on the relative momentum angle or explicitly on φ x . Therefore, we can shift all four integration momentum angles by φ x and eliminate φ x everywhere except for the explicit terms. After doing the φ x integral the cross-terms in Eq. (A.7) vanish, while terms with cos 2 φ x and sin 2 φ x add up to 1 2 cos(φ ∆x − φ p ). Since the integral is symmetric in p, k, p and k , except for p ⊥ and φ p in cos(φ ∆x − φ p ), we will symmetrize over all four momentum variables and divide the integral by 4. Introducing a shorthand notation we arrive at our final expression The multidimensional integral present in the linearized form of Eq. (A.9) has been evaluated numerically using Monte Carlo with importance sampling for discrete values of r and τ . For the explicit phase-space parametrization of the integral see ref. [51]. Mean values and stochastic errors presented in the figures and tables throughout this article have been calculated using the jackknife resampling method. Discretization errors in the r − τ plane with a N r × N τ = 20 × 20 grid were checked to be of percent size and negligible in comparison to the statistical uncertainties.
In order to verify the results of both analytical predictions and numerical simulations, a number of crosschecks have been preformed. The validity of the derivation and implementation of the linerized expression of D n , Eq. (A.9), has been checked against a nonlinear implementation, i.e. evaluating Eq. (2.11) directly. However, the nonlinear method has larger statistical uncertainties, which become worse for very small values of . Therefore all reported results are obtained with the linearized equation Eq. (A.9). We have also implemented RTA kinetic theory Eq. (3.13) in the same setup. As both kinetic theories share the same free-streaming background evolution, the reproducability ITA results in Eq. (3.15) was used as additional check.

B Results of the parameter space scan
To numerically test the scaling relations predicted in Section 3.1 we performed a systematic variation of all model parameters: starting time τ 0 , Gaussian width of the density profile R 0 , the normalization of the distribution A, anisotropy parameter ξ and the coupling constant λ. We tabulate the values of the screening mass m 2 g at initial time at the origin (r 0 = 0.01), scaling variablem 2 g and the classical and Bose-enhanced contributions to the elliptic flow. For the linearized approach, we can completely scale the value of the eccentricity ε n . The results for the CGC-like momentum distribution, Eq. (2.20), are given in Table 1. In addition, we performed simulations with deformed thermal initial conditions, Eq. (2.26) and the results are given in Table 2. Finally, for completeness in Table 3 we record the results of simulations for which we do not expect scaling, because τ 0 /R 0 1 assumption is violated. These results were used to produce Figs. 1 and 2.

C Initial conditions in nuclear collisions
In order to determine realistic values of the scaling variables corresponding to a physical collision systems, we will re-use the tabulated values of eccentricity ε 2 , RMS entropy radius r 2 and entropy density dS/dy/(π r 2 ) (in arbitrary units) from ref. [47]. These initial conditions for PbPb, OO, pPb and pp collision systems were generated using the TRENTo initial state model [48,49]. The entropy normalization is not specified, therefore we will use the total entropy per rapidity dS/dy = 11335 extracted from data for √ s NN = 2.76 TeV PbPb 0-10% collisions [52]. It is known experimentally that the particle multiplicity (∝ entropy) in nucleus-nucleus collisions scales with the collision energy according to ∝ ( √ s) 0.31 law [53]. Therefore we use a (5.02/2.76) 0.31 factor to increase the entropy. The next step is to convert this final state entropy to the initial state energy density. In homogeneous boost-invariant systems the early time non-equilibrium evolution of energy density can be well described by hydrodynamic attractor curves. We will use the following formula derived in ref.
[54] to relate entropy density per rapidity (sτ ) final to initial energy density per rapidity (eτ ) 0 (eτ ) 0 = (sτ ) Here C ∞ ≈ 0.9 is the property of the hydrodynamic attractor in QCD EKT, η/s is the specific shear viscosity and ν eff ≈ 40 is the effective number of degrees of freedom in the high temperature equilibrium QGP phase. Note that the QCD EKT attractor was obtained for an initial state with only gluons present, so we can identify (eτ ) 0 as the energy density of the gluonic degrees of freedom. We assume that the averaged initial energy density profile in each centrality class is described by a Gaussian (eτ ) 0 exp(−|x ⊥ | 2 /R 2 0 ), where R 0 is 2/3 times the RMS of the entropy profile in ref. [47]. Then the energy density at the origin (eτ ) 0 is obtained using Eq. (C.1) with η/s = 0.2.