Event-by-event fluctuations in the medium-induced jet evolution

We develop the event-by-event picture of the gluon distribution produced via medium-induced gluon branching by an energetic jet which propagates through a dense QCD medium. A typical event is characterized by the production of a large number of soft gluons which propagate at large angles with respect to the jet axis and which collectively carry a substantial amount of energy. By explicitly computing 2-gluon correlations, we demonstrate the existence of large event-by-event fluctuations, which reflect the stochastic nature of the branching process. For the two quantities that we have investigated — the energy loss at large angles and the soft gluon multiplicity —, the dispersion is parametrically as large as the respective expectation value. We identify interesting scaling laws, which suggest that the multiplicity distribution should exhibit KNO (Koba-Nielsen-Olesen) scaling. A similar scaling is known to hold for a jet branching in the vacuum, but the medium-induced distribution is found to be considerably broader. We predict that event-by-event measurements of the di-jet asymmetry in Pb+Pb collisions at the LHC should observe large fluctuations in the number of soft hadrons propagating at large angles and also in the total energy carried by these hadrons.


Introduction
A remarkable phenomenon discovered and, by now, abundantly studied in Pb+Pb collisions at the LHC is the di-jet asymmetry -a strong energy imbalance between two jets which propagate nearly back-to-back in the plane transverse to the collision axis [1][2][3][4][5][6][7][8][9]. Imbalanced di-jets are also observed in proton-proton collisions, where they are mostly associated with 3-jet events, but the respective events in (central) Pb+Pb collisions look quite different. First, the fraction of events showing a strong energy asymmetry is considerably higher. Second, the pattern of the asymmetry, i.e. the distribution of the 'missing' energy in the underlying event, is very different: in heavy ion collisions, the energy difference between JHEP05(2016)008 the leading and the subleading jets is not carried by a 3rd jet anymore, but rather by many soft hadrons, with energies p T ≤ 2 GeV, which propagate in the hemisphere of the subleading jet at large angles w.r.t. the latter. Accordingly, the 'missing' energy is only slowly recovered when increasing the angular opening ∆R of the subleading jet; for instance, for the relatively large value ∆R = 0.8, there is still a substantial energy imbalance, ∆E 20 ÷ 30 GeV [2,8].
It seems natural to associate the di-jet asymmetry in Pb+Pb collisions with the interactions between the partons in the jets and the surrounding QCD medium -the quark-gluon plasma -which is expected to be created in the intermediate stages of a ultrarelativistic nucleus-nucleus collision. Such interactions (e.g. elastic collisions) can easily deviate the soft components of a jet towards large angles and an asymmetry is generated whenever one of the two jets crosses the medium along a longer distance than the other [10]. However, the main question remains: how is that possible that a substantial fraction of the energy of a jet is carried by its soft constituents ? As well known, bremsstrahlung in QCD cannot lead to an efficient redistribution of the energy towards soft particles: while this mechanism favors the emission of soft gluons, the total energy carried by these gluons remains negligible. This is so since each radiated gluon carries only a tiny fraction x 1 of the energy of its parent parton.
However, as pointed out in ref. [11], the evolution of a jet in a dense medium is dramatically different from that that would be generated by bremsstrahlung in the vacuum: at sufficiently low energies, the BDMPSZ rate for medium-induced gluon branching favors multiple branching which is quasi-democratic. The BDMPSZ rate (from Baier, Dokshitzer, Mueller, Peigné, Schiff, and Zakharov) [12][13][14][15][16] refers to gluon emissions which are triggered by multiple soft scattering between the gluon which is about to be emitted and the constituents of the medium (see also [17][18][19][20][21] for related developments). Unlike the rate for bremsstrahlung, the BDMPSZ rate depends also upon the energy of the parent parton (and not only upon the energy fraction x of the emitted gluon). Namely, this rate becomes large when the emitter is sufficiently soft (irrespective of the value of x) and thus favors the development of democratic cascades via the successive branchings of soft gluons. The democratic branching is a very efficient mechanism for transmitting the energy from the leading particle (the energetic parton which has initiated the jet) to a myriad of soft gluons, which are eventually deviated to large angles by rescattering in the medium.
As demonstrated in refs. [22][23][24], the medium-induced evolution via multiple branching is a classical stochastic process, which is Markovian: 1 after each splitting, the two daughter gluons evolve independently from each other because their color coherence is rapidly washed out by rescattering in the medium. So far, only the mean field aspects of this stochastic evolution, as encoded in the average gluon spectrum D(x) = x(dN /dx) and the associated transport equation, have been explicitly studied. For the branching dynamics alone, one has obtained exact solutions for D(x) [11,37,38], which demonstrate the physics of democratic 1 This property was implicitly assumed in earlier works aiming at constructing a kinetic theory for in-medium jet evolution [25][26][27][28][29]. The study of color decoherence in a dense partonic medium has been pioneered in refs. [30][31][32] and the precise conditions for the factorization of successive emissions have been clarified in refs. [22,24]. See also refs. [33][34][35] for related studies and the recent review paper [36].

JHEP05(2016)008
branchings and the associated phenomenon of wave turbulence (the energy flux down the democratic cascades is uniform in x). Using more general equations which include the elastic collisions with the medium constituents, one has also studied the interplay between branchings and elastic collisions, and its consequences for the angular structure of the cascade [23,34,[39][40][41] and for the thermalization of the soft branching products [42].
In this paper, we shall study the event-by-event fluctuations introduced by the branching process in the medium-induced jet evolution. Specifically, we shall compute the dispersion in the energy lost by the jet at large angles and also in the number of gluons N (x 0 ) having an energy fraction x larger than some infrared cutoff x 0 1. Here, x ≡ ω/E, with E the energy of the leading particle and ω ≤ E the energy of a parton from the jet. We shall consider a pure branching process, without elastic collisions, so our results strictly apply only for sufficiently high energies ω T , with T the typical energy of the medium constituents. (The medium is assumed to be a weakly-coupled quark-gluon plasma in thermal equilibrium at temperature T ; for the present purposes, it is fully characterized by a transport coefficient for transverse momentum diffusion, known as the 'jet quenching parameter'q.) This is not a serious limitation so long as we are interested only in the overall energy transmitted by the jet to the medium, and not also in its detailed distribution in space and time [42]. For this branching process, we shall obtain an exact result, eq. (3.9), for the gluon pair density D (2) (x, x ) ≡ xx (dN/dxdx ). Using this result, we shall compute the aforementioned dispersions.
Our main conclusion is that fluctuations are large: for both the energy loss at large angles and the gluon multiplicity N (x 0 ), the dispersion is parametrically as large as the respective average quantity. Such large fluctuations should be easy to observe in eventby-event studies of the di-jet asymmetry at the LHC. In particular, we predict that the multiplicity fluctuations should exhibit Koba-Nielsen-Olesen (KNO) scaling [43]. A similar scaling is known to hold for a jet branching in the vacuum [44], but the medium-induced gluon distribution is found to be considerably wider (see section 4 for details). The importance of fluctuations for the phenomenology of di-jet asymmetry is also demonstrated by a very recent numerical study [45], which appeared when our paper was about to be finalized. We shall further refer to this study in our concluding section.
The physical picture emerging from our analysis can be summarized as follows. In a typical event, the leading particle evolves by emitting a number of order one of primary gluons with energy ω ∼ ω br (L) together with a large number ∼ [ω br (L)/ω] 1/2 of considerably softer gluons, with ω ω br (L). Here, ω br (L) ≡ᾱ sq L 2 , withᾱ s = α s N c /π and L the distance travelled by the jet through the medium, is the characteristic energy scale for medium-induced multiple branching. (We implicitly assume here that ω br (L) E, as this is the typical situation for high energy jets at the LHC; see the discussion at the end of section 2.1.) The primary gluons with ω ω br (L) develop mini-jets via successive democratic branchings and thus transmit their whole energy to softer quanta with ω ∼ T that are expected to thermalize via collisions in the medium [42]. Harder emissions with ω ω br (L) occur only rarely, with a small probability ∼ [ω br (L)/ω] 1/2 , and moreover they do not contribute to the energy lost by the jet as a whole, since hard gluons propagate at small angles. The energy loss at large angles is rather controlled by the hardest typical JHEP05(2016)008 emissions, those with energies ω ∼ ω br (L). As aforementioned, the number of such emissions is of order one and they occur independently from each other; accordingly, both the average energy loss by the jet and its dispersion are of order ω br (L) (see eq. (3.15)).
Unlike the primary gluons (the direct emissions by the leading particle), which are quasi independent from each other, the gluons from the subsequent generations, as produced via democratic branchings, can be mutually correlated due to the fact that they have common ancestors. This correlation is particularly important for the sufficiently soft gluons, with ω ω br (L), which belong to a same mini-jet. It is responsible for the strong fluctuations in the gluon multiplicity at small x, leading to the KNO scaling alluded to above (see eq. (4.10) and the subsequent discussion).
Another interesting conclusion of our analysis is that the gluon distribution at small x is not affected by multiple branching: both the gluon spectrum D(x) and the gluon pair density D (2) (x, x ) are formally the same, when x, x 1, as the respective predictions of the lowest-order perturbative expansion in the number of branchings. For the gluon spectrum, this property was already known [11,25]: this is the statement that the low-x limit of the BDMPSZ spectrum, namely the power-law spectrum D(x) ∝ 1/ √ x, is a fixed point of the branching dynamics. The present analysis shows that a similar property holds for the 2-gluon correlation D (2) (x, x ) (see eq. (3.10) and the related discussion). This might look surprising for the soft gluons with ω ω br (L), for which multiple branching is known to be important, but it is most surely a consequence of wave turbulence: the rate at which the soft gluons are produced via the splitting of harder gluons is precisely equal to the rate at which they disappear by decaying into even softer gluons. This fine balance between 'gain' and 'loss' ensures that the energy flux is independent of x, which is the hallmark of wave turbulence [46,47].
While not manifest in the gluon distribution at small x, the effects of multiple branchings are visible in other aspects of this distribution, like the structure of the leading particle peak near x = 1, and also in global properties, like the phenomenon of wave turbulence and its consequences for the energy loss by the jet at large angles. Such phenomena are quite intricate and often elusive, hence the importance of disposing of exact analytic solutions in order to unambiguously demonstrate them. Analytic approximations and parametric estimates can be helpful for developing a general picture, but they cannot capture subtle aspects like the broadening of the leading particle peak via multiple emissions, or the existence of a turbulent energy flux, uniform in x. Numerical solutions too may alter some aspects of the dynamics, like the power-law spectrum at small x, due to artifacts like infrared cutoffs. This being said, and once a coherent physical picture has been established on the basis of exact results, the numerical methods offer a powerful tool for extending these solutions to more realistic situations and, especially, for systematic studies of the phenomenology. In particular, the use of Monte-Carlo implementations of the in-medium jet evolution like MARTINI [29] and JEWEL [48] should allow for general studies of the stochastic aspects that we shall identify in this paper and of their phenomenological implications. This paper is organized as follows. In section 2 we succinctly describe the physical picture of the medium-induced jet evolution and its mathematical formulation as a Marko-JHEP05(2016)008 vian process. In particular, we present the transport equation (2.14) obeyed by the gluon pair density D (2) (x, x ). More details on the formalism are deferred to appendix A. In section 3 we discuss the energy loss at large angles, operationally defined as the total energy transmitted, via successive branchings, to the very soft gluons with x → 0. Section 3.1 is devoted to the mean field picture, that is, the gluon spectrum D(x) and the average energy loss. Most of the results presented there were already known, but our physical discussion is more furnished, in line with our general purposes. In section 3.2, we present our main new results, which are both exact (within our theoretical framework): eq. (3.9) for the gluon pair density D (2) (x, x ) and eq. (3.13) for the variance in the energy loss at large angles. The physical interpretation of these results is discussed at length, in section 3.2 and the dedicated section 3.3. Details on the calculations are presented in appendices B and C. In section 4, we discuss the gluon number distribution, for gluons with energy fraction x ≥ x 0 . In section 4.1 we compute the average multiplicity, while in section 4.2 we present and discuss our results for the second factorial moment N (N − 1) and for the variance. The respective calculations are quite tedious (the details are deferred to appendix D), but in section 4.2 we anticipate the final results via physical considerations, which also shed light on their remarkable scaling properties. Section 5 presents a brief summary of our results together with possible implications for the phenomenology.
2 Medium-induced jet evolution: the general picture As explained in the Introduction, we consider the jet generated via medium-induced gluon branching by a leading particle (LP) with initial energy E which propagates along a distance L through a weakly coupled quark-gluon plasma with temperature T E. Some typical values, as inspired by the phenomenology at the LHC, are E = 100 GeV and T = 0.5 GeV.
We would like to study the event-by-event distribution of the energy lost by the jet at large angles with respect to the jet axis, say at polar angles θ > θ 0 with θ 0 ∼ O(1).
To that aim, one needs the distribution of the radiated gluons in energy (ω) and polar angle (θ), or, equivalently, in ω and transverse momentum k ⊥ (recall the trigonometric relation sin θ = k ⊥ /ω). For simplicity though, and in order to allow for analytic studies, we shall explicitly consider only the distribution in ω, as obtained after integrating out k ⊥ . Then the angular distribution will be approximately reconstructed by using the fact that, for medium-induced emissions, the k ⊥ -distribution is strongly peaked at 2 k 2 ⊥ = Q 2 s ≡qL. Indeed, this is the typical transverse momentum acquired via medium rescattering by a gluon propagating through the medium along a distance L. For the phenomenologically relevant valuesq = 1 GeV 2 /fm and L = 4 fm, one finds Q s 2 GeV, meaning that gluons with energies ω ≤ 2 GeV will propagate at angles θ 1.
2 This estimate for k 2 ⊥ ceases to be valid for the very soft gluons with energies ω ω br (L), whose lifetimes t br (ω) are much shorter than L (see eq. (2.4)). For such gluons, k 2 ⊥ qt br (ω) is smaller than Q 2 s and also ω-dependent [34,39,40]. However, this brings no serious complication for our present arguments, since the energy loss in eq. (2.1) is quasi-independent of x0, hence of ω0, whenever ω0 ω br (L); see figure 3 and the associated discussion.

JHEP05(2016)008
The above considerations motivate the following estimate for the energy which propagates at angles θ > θ 0 at time t, with 0 ≤ t ≤ L: Here, dN /dω is the energy distribution of the gluons from the jet in a given event, 3 x ≡ ω/E is the energy fraction of a gluon, x 0 ≡ ω 0 /E, and the infrared cutoff ω 0 is the energy corresponding (via the simple argument above) to the limiting angle θ 0 : sin θ 0 = Q s /ω 0 . The r.h.s. of eq. (2.1) is recognized as the difference between the initial energy E of the LP and the energy carried by the components of the jet whose energies are larger than ω 0 (and hence propagate at angles smaller than θ 0 ). For the jets produced in actual nucleus-nucleus collisions, this quantity E(t, x 0 ) can significantly vary from one event to another, due to various reasons. First, the medium properties (like the value ofq) and its geometry (size and shape) can be different in different events -albeit in practice one can diminish the importance of such fluctuations by choosing events in a same centrality class. Second, even for a given collision geometry, there can be large fluctuations in the distance L travelled by the jet through the medium, because the LP can be produced at any point within the interaction volume and with any initial orientation. Finally, the evolution of the jet via successive medium-induced branchings is a stochastic process, which is characterized by strong fluctuations, as we shall see.
In what follows, we shall only address the last source of fluctuations, that is, we shall assume fixed values for the physical parametersq, L and E, and study the distribution of E(x 0 , t) generated by the stochastic branching process. To that aim, we shall rely on the Markovian description for this process, as recently developed in refs. [11,23,[36][37][38]. In this section we shall succinctly review this description, starting with the underlying physical picture. Besides motivating the various approximations, this discussion will be also useful for the physical interpretation of our new results.

Democratic branchings and wave turbulence
Consider a generic gluon from the medium-induced cascade, whose energy ω satisfies T ω ≤ E. This gluon has a probability to radiate a soft gluon with energy ω ≥ zω during a time interval ∆t (z 1 is the splitting fraction). In eq. (2.2),ᾱ s ≡ α s N c /π,q is the jet quenching parameter, and t f (zω) = zω/q is the quantum formation time for the emission of a gluon with energy zω. So, eq. (2.2) can be understood as follows: ∆P is the product of the probabilityᾱ s for one gluon emission JHEP05(2016)008 times the number of times that this emission can fit within a given time interval ∆t. This probability becomes of order one, meaning that a branching is certain to occur, after a time This time is relatively short when z 1, meaning that soft gluons are abundantly produced. Clearly, the original gluon survives to such soft emissions, in the sense that it can be unambiguously distinguished from its radiation products, due to its higher energy. But after the larger time interval the gluon ω is bound to undergo a 'quasi-democratic' (z ∼ O(1)) branching and then it disappears -it gets replaced by a pair of daughter gluons whose energies are comparable to each other. In turn, these daughter gluons disappear (via democratic splittings) even faster, since their energies are lower. This democratic branching process repeats itself until the original energy ω gets transmitted to a multitude of soft gluons with comparable energies. This happens over a time interval of the order of t br (ω). Eventually, the dynamics changes when the energies of the soft descendants become of the order of the medium scale T : such very soft gluons can efficiently thermalize via elastic collisions off the medium constituents and thus they transmit their whole energy to the medium [42]. A fundamental property of a democratic cascade, known as 'wave turbulence', is the fact that the energy flux generated via democratic branchings -the rate at which the energy flows across a given bin in ω -is independent of ω, for ω E [11,37]. Through successive branchings, the energy flows from one parton generation to the next one, without accumulating at any intermediate value of ω. In the absence of thermalization effects, the whole energy would eventually accumulate into a condensate at ω = 0 [11]. In reality, the energy is eventually transmitted to the medium constituents, via elastic collisions leading to the thermalization of the soft branching products with ω T [42]. But the rate at which this happens is still controlled by the branching process at ω > T , that is, by the turbulent flux alluded to above, which is independent of ω and hence can be formally computed at ω = 0.
These considerations suggest that, in order to evaluate the energy loss in eq. (2.1), one can restrict oneself to a pure branching process (no elastic collisions) and take the limit x 0 → 0. This is consistent with the fact that the energy which is physically lost towards the medium is the same as that which would accumulate at x = 0 in the absence of thermalization. Hence, one can replace eq. (2.1) by with dN /dx determined by the pure branching process discussed at length in refs. [11,23,[36][37][38].
We conclude this subsection with some remarks of relevance for the phenomenology of high energy jets (E ≥ 100 GeV) at the LHC. In discussing democratic branchings so JHEP05(2016)008 far, we have implicitly assumed that the branching time (2.4) is shorter than the medium size L. This in turn implies an upper bound on the energy ω of the parent gluon: ω ω br (L) ≡ᾱ 2 sq L 2 . This condition however is not satisfied by the LP in the experimental set-up at the LHC: its energy E ≥ 100 GeV is considerably larger than the typical values expected for ω br (L). (For instance, withᾱ s = 0.3,q = 1 GeV 2 /fm, and L = 4 fm, one finds ω br 8 GeV.) This means that the LP cannot undergo (medium-induced) democratic splittings and hence it emerges out of the medium, as the core of the surviving jet. A similar argument applies to the relatively hard emissions 4 with energies ω ω br (L), which are rare events: such hard gluons survive the medium and propagate at small angles, hence they appear in the final state as components of the conventionally defined jet.
On the other hand, the LP can emit abundantly, i.e. with probability of order one, relatively soft gluons with ω ω br (L). Such soft primary emissions determine the structure of a typical event. Each of these primary gluons will generate its own cascade, or 'mini-jet', via successive democratic branchings and thus eventually transmit its initial energy to the medium. The 'hardest' among these soft gluons, those with energies of order ω br (L), will play an essential role in what follows: they control both the energy loss at large angles and its fluctuations, as we shall see.

Transport equations for the gluon distribution
In this section, we shall relate the average energy loss E(t) and its variance E 2 (t) − E(t) 2 to the one-and two-point functions of the energy density x(dN /dx), for which we shall then establish evolution equations. For the average energy loss, this relation is quite obvious: where ε ≡ E/E is the fraction of the total energy which accumulates at x = 0 and D(x, t) is the average energy density, or gluon spectrum, as defined in eq. (2.9) below. For the variance, we have As we shortly explain, the quantity X 2 (t) can be evaluated as where D (2) (x, x , t), with x + x ≤ 1, is the average density of pairs of gluons multiplied by xx : The energy of a medium-induced emission is limited by the condition that the formation time t f (ω) be smaller than L. This implies ω ≤ ωc(L) ≡qL 2 [12,14]. However, the relatively hard gluons with ω br (L) < ω < ωc(L) cannot undergo democratic branchings and hence do not contribute to the energy loss at large angles.

JHEP05(2016)008
To render these definitions more explicit, let us consider a more detailed description of the statistical ensemble of events which is created by the stochastic branching process [23]. The central ingredient is the probability density P n (x 1 , x 2 , · · · , x n |t) for having a state with n gluons with energy fractions x i (i = 1, . . . , n), at time t. At t = 0, we have just the LP, hence P n (t = 0) = δ n1 δ(x 1 − 1). The expectation value of an arbitrary observable is computed as where O n ≡ O(x 1 , x 2 , · · · , x n ) denotes the value of O in a particular event with n gluons.
In particular, the quantities introduced in eq. (2. Note the distinction i = j in the definition of D (2) (x, x , t), which shows that this quantity is indeed the density of pairs of gluons, each pair being counted twice (since the pairs (x, x ) and (x , x) are separately counted). Using eq. (2.11), it is easy to check the previous formula (2.8).
At this point, it is important to observe that this probabilistic description requires an infrared cutoff (say, a lower limit on x), playing the role of an energy resolution scale, below which gluons cannot be resolved anymore. Indeed, the branching dynamics produces an infinite number of arbitrarily soft gluons and the 'state with exactly n gluons' is not well defined without such a cutoff. Any explicit construction of such a state, say via Monte-Carlo simulations, must involve an infrared cutoff on x, to be viewed as a part of the 'state' definition. On the other hand, physical quantities which are sufficiently inclusive, such as D(x, t) and D (2) (x, x , t), are insensitive to the unobserved, soft, gluons, hence they must be independent of this cutoff. And indeed, one can show that the cutoff dependence cancels out in the equations obeyed by these quantities. So long as one is solely interested in such quantities, one can formally proceed without any infrared cutoff (see the discussion in appendix A).
The probability densities P n (x 1 , x 2 , · · · , x n |t) obey 'master equations', that is, coupled rate equations whose general structure is rather standard, since common to all the Markovian branching processes. (These equations are presented in appendix A for convenience; see also [23].) What is specific for the problem at hand, is the expression of the splitting rate, here given by the BDMPSZ mechanism [12,14]. Namely, the differential probability per unit time and per unit z for the splitting of the parent gluon ω = xE into the pair of daughter gluons zω and (1 − z)ω can be conveniently written as with t br (ω) the branching time from eq. (2.4). Using the master equations within eq. (2.10), it is straightforward to derive transport equations for the observables (see appendix A for JHEP05(2016)008 details). The corresponding equation for the gluon spectrum D(x, t) reads 5 where we introduced the dimensionless time variable τ ≡ t/t br (E), that is, the physical time in units of the branching time t br (E) of the LP. The first term in the r.h.s. describes the gain in the number of gluons at x due to emissions from gluons with x = x/z > x, whereas the second term describes the loss via the decay into softer gluons. Taken separately, the gain term and the loss term have endpoint singularities at z = 1 (corresponding to the emission of very soft gluons), but these singularities mutually cancel and the overall equation is well defined without any infrared cutoff, as anticipated. The corresponding equation for the two-point function will be derived in appendix A, and reads 14) The two terms in the r.h.s. of the first line are easy to understand by comparison with eq. (2.13): they describe the separate evolutions in the numbers of gluons with energy fractions x and x , respectively. The 'source' term in the second line describes the simultaneous creation of a pair of gluons with energy fractions x and x via the branching of a parent gluon with energy fraction x + x (with x + x ≤ 1 of course). This terms introduces correlations in the evolution of the pair of gluons and is responsible for the fluctuations to be analyzed in the next sections. The initial conditions for the above equations read D(x, 0) = δ(x − 1) and D (2) (x, x , 0) = 0. The respective solution D(x, τ ) to eq. (2.13) acts as a Green's function for eq. (2.14): the solution to the latter can be written as where S(x 1 , x 2 , τ ) is a compact notation for the source term in the r.h.s. of eq. (2.14). The physical interpretation of eq. (2.15) is quite clear: at some intermediate time τ , a gluon with energy fraction x 1 + x 2 splits into two daughter gluons with energy fractions x 1 and respectively x 2 , whose subsequent evolutions generate the mini-jets within which one measures the two final gluons x and x . Note that the parent gluon with energy x 1 + x 2 is the last common ancestor (LCA) of the two measured gluons x and x : after this ancestor has split, the subsequent branchings leading to x and x are completely disconnected from each other (see figure 1 for an illustration). 1, have a last common ancestor (LCA) with energy fraction x 1 + x 2 , whose splitting vertex is indicated by a blob. Left: the LCA is a primary gluon with x 1 + x 2 1 and the two measured gluons belong to a same democratic cascade. Right: the LCA is the LP itself (x 1 + x 2 1), which emits a primary gluon with x 2 1; in this case the two measured gluons belong to different mini-jets.

Event-by-event fluctuations in the energy loss at large angles
We are now prepared for an explicit calculation of the two quantities of primary interest for us here: the average energy lost by the jet at large angles E(t) and its variance To obtain exact analytic results, we shall consider a slightly simplified version of the stochastic process introduced above, as obtained by replacing the kernel K(z) in eq. (2.12) with K 0 (z) ≡ 1/[z(1 − z)] 3/2 . This replacement is pretty harmless: the simplified kernel preserves the poles of the exact kernel at z = 0 and z = 1 and the respective residues, so it generates a very similar evolution. And indeed, the numerical solutions to eq. (2.13) using the exact kernel [37,38] are very similar, even quantitatively, to the exact analytic solution presented in ref. [11] for the simplified kernel K 0 (z).

The gluon spectrum and the average energy loss
The solution to eq. (2.13) with K(z) → K 0 (z) and the initial condition D(x, 0) = δ(x − 1) reads [11] This solution is illustrated in figure 2. Its physical interpretation is important for what follows, so we shall now elaborate on it. To that aim, let us first recall that where ω br (t) =ᾱ 2 sq t 2 is the energy of the hardest emissions that can occur with a probability of O(1) during a time t. The most interesting regime for us here is the small-τ regime τ 1, which corresponds to the typical experimental situation at the LHC: a very energetic LP JHEP05(2016)008 with E ω br (L) or, equivalently, a relatively small medium size L t br (E). We recall that L is the total distance travelled by the LP through the medium and hence the maximal value of t.
At small times, such that πτ 2 1, eq. (3.1) exhibits a pronounced peak just below x = 1, which describes the leading particle. The position x p (τ ) of the peak (the local maximum of the function D(x, τ )) represents the most likely value of the energy of the LP at time τ . Hence, the difference 1 − x p is the most likely value for its energy loss, that is, the energy loss by the LP in a typical 6 event. One easily finds 1 − x p (2π/3)τ 2 , or, in physical units, ∆E ≡ E(1 − x p ) ∼ ω br (t), in agreement with the discussion in section 2.1: during a time t, the LP radiates with probability of O(1) gluons with energies ω ω br (t). Its energy loss is dominated by the hardest among these emissions. The broadening of the LP peak, say, as measured by its width δx at half the height of the peak, is itself of order πτ 2 . Clearly, this width δx is a measure of the dispersion in the energy distribution of the LP. The fact that the dispersion is comparable with the mean energy loss, δx ∼ 1−x p ∼ τ 2 , is an interesting observation, whose interpretation will be discussed in section 3.3.
Consider also the tail of the gluon spectrum at small x, which describes the soft radiation. In view of the previous discussion, one might expect an accumulation of gluons in the small-x bins at x πτ 2 , but this is actually not the case: for x 1, eq. (3.1) exhibits the same power-law tail D(x, τ ) ∝ 1/ √ x as the BDMPSZ spectrum generated via JHEP05(2016)008 a single emission by the LP: with the second approximation valid at small times πτ 2 1. In other terms, the spectrum at small x is not modified by multiple branching. Mathematically, this is possible because the power-law spectrum in eq. (3.3) represents a fixed point 7 for the rate equation (2.13) [11,25]: the 'gain' and 'loss' terms precisely cancel each other. This condition is necessary to ensure that the energy flux associated with multiple branching is independent of x ('wave turbulence'); see [37] for a calculation of this flux.
At later times πτ 2 1, the LP peak disappears and the function D(x, τ ) has support at x 1/πτ 2 . With increasing τ , both the support in x and the strength of the spectrum at any point x within this support are rapidly decreasing, meaning that the energy flows outside the spectrum. To characterize this flow, we now compute the average energy fraction carried by the whole parton cascade at time τ . Using eqs. (2.6) and (3.1), one easily finds [11] where in performing the integral it was convenient to change variable as The integral is controlled by u 1/πτ 2 , meaning by x 1/(1 + πτ 2 ). At small times πτ 2 1, this yields x 1 − πτ 2 , showing that, not surprisingly, most of the jet energy is carried by the LP. What is more interesting though, is that the missing energy at small times, namely 5) or, in physical units (cf. eq. (3.2)), is comparable with the energy radiated by the LP in a typical event (recall the discussion after eq. (3.1)). This is consistent with the above observation that there is no energy accumulation in gluon spectrum at x > 0: via successive democratic branchings, the primary gluons with x τ 2 transmit their energy to arbitrarily soft quanta, which formally accumulate at x = 0. (Physically, they thermalize and thus transmit their whole energy to the medium [42].) Equations like eq. (3.3) and eq. (3.4) also confirm the fact that the lifetime of the jet as a whole is of order t br (E), i.e. it is comparable with the time it takes the LP to undergo a first democratic branching. It is also interesting to study the energy carried by the gluons in the spectrum having an energy fraction larger than a given value x 0 1. Via the relation sin  propagates inside a cone with opening angle θ 0 . (Recall the discussion around eq. (2.1).) The corresponding generalization of eq. (3.4) reads where u 0 ≡ x 0 /(1 − x 0 ) and we have introduced the error function The last, approximate, equality in eq. (3.7) holds when τ √ x 0 1. In figure 3, we exhibit X(τ, x 0 ) as a function of 1/x 0 for several interesting values of τ ≤ 1. As manifest in this figure, X(τ, x 0 ) is essentially flat at 1/x 0 1. In other terms, the total energy contained inside the jet increases only slowly when increasing the jet angle θ 0 ∝ 1/x 0 . This is easy to understand in the present context: when x 0 1, the energy contained in the soft tail of the spectrum at 0 < x < x 0 is very small and truly negligible not only compared to the energy carried by the harder gluons with x x 0 , but also compared to the 'missing' energy ε that has accumulated at x = 0. (In figure 3, this missing energy is represented by the 'offset' 1 − X(τ, x 0 ) of the various curves at 1/x 0 → ∞.) This discussion reenforces our previous argument, around eq. (2.5), that the energy loss at large angles can be accurately computed by letting x 0 → 0 in eq. (2.1). In principle, one needs a non-zero cutoff x 0 ∝ 1/θ 0 , to control the angular distribution of the radiation. However, the values of x 0 corresponding to large angles θ 0 ∼ 1 are quite small, JHEP05(2016)008 x 0 1, and for them the x 0 -dependence of the energy distribution is truly negligible, as shown by figure 3. This ultimately reflects the fact that the energy loss at large angles is controlled by the turbulent energy flow, which is independent of x.

The gluon pair density and the variance in the energy loss
We now turn to our main new calculation, that of the variance σ 2 ε (τ ) in the energy loss at large angles. As explained in section 2.2, this involves the gluon pair density D (2) (x, x , τ ), for which an exact but formal solution has been written in eq. (2.15). Remarkably, it turns out that for the simplified kernel K 0 (z) (i.e., for the gluon spectrum D(x, τ ) in eq. (3.1)), the integrals in eq. (2.15) can be computed exactly, with the following result (see appendix B for details): To get more insight into this result, consider the particular case where πτ 2 1 (that is, the LP has no time to undergo a democratic branching) and x, x 1. In this limit, eq. (3.9) reduces to where we have also used eq. (3.3). This result is truly remarkable in several respects. First, eq. (3.10) coincides with the respective prediction of lowest-order perturbation theory, 8 that is, the 3 graphs in figure 4 which describe the emission of only 2 soft gluons. (Indeed, one needs at least 2 branchings in order to produce 2 soft gluons in the final state.) These 3 graphs can be easily evaluated and one finds that they give identical contributions; in particular, the overall factor 1/2 comes in each of them from integrating over the branching time τ of the last common ancestor (LCA): A priori, this perturbative estimate is expected to be correct for sufficiently hard emissions, x, x τ 2 , but the result in eq. (3.3) also applies to very soft gluons, with x, x τ 2 , for which the effects of multiple branchings are known to be importantactually, non-perturbative. This observation extends to the gluon pair density an important JHEP05(2016)008 property that we noted in section 3.1 in relation with the gluon spectrum (recall the discussion after eq. (3.3)): the fact that the gluon distribution produced by the branching process at small x is formally insensitive to multiple branching, since controlled by the turbulent fixed point. Of course, the effects of multiple branching become visible when looking at the structure of the leading-particle peak (see e.g. eq. (3.9) for x + x 1) and also at global properties, like the total energy enclosed in the spectrum. But, at least for small times τ 1 (so long as the LP still exists), they are not manifest in the gluon distribution at small x.
The physical interpretation of the result in eq. (3.3) for generic values x, x 1 can be inferred from the explicit calculation of the integrals in eq. (2.15), in appendix B. This involves 3 types of processes which are topologically similar to those depicted in figure 4 -the latter must be simply dressed by multiple branchings.
In one of these processes, the two measured gluons x and x belong to a same mini-jet (see figure 1 left). That is, their LCA, which has energy fraction x 1 +x 2 in eq. (2.15), is itself soft (x 1 + x 2 1), but such that x 1 − x ∼ τ 2 and similarly x 2 − x ∼ τ 2 . These conditions ensure that the gluon x is found with probability of O(1) in the cascade generated by the gluon x 1 during a time interval τ − τ ∼ τ , and similarly for the pair of gluons x and x 2 . If the measured gluon, say x, is relatively hard, such that x τ 2 (with x 1 though), then x 1 x, meaning that the energy fractions x and x 1 refer to a same gluon, which emits some soft radiation in passing from x 1 to x. If on the other hand x τ 2 , then x, so the measured gluon x is one of the soft gluons in the cascade generated by x 1 via democratic branchings. A similar discussion applies to x and x 2 .
In the 2 other processes, the LCA is the LP (x 1 + x 2 1), so the two measured gluons belong to different mini-jets. E.g. the LP emits a primary gluon x 2 with x 2 − x ∼ τ 2 , which then generates a democratic cascade which includes the final gluon x . The other final gluon, with energy x, is found in the mini-jet produced by another primary gluon, which is later radiated by the LP (see figure 1 right).
) is a measure of the correlations in the gluon distribution. We shall later argue, in section 3.3 and in section 4, that the successive emissions of primary gluons are quasi-independent from each other. This in turn implies that the net correlation D (2) is to be associated with the first process described above, cf. figure 1 left, where the two measured gluons belong to a same mini-jet.
Using eq. (3.9), it is straightforward to compute the first piece in the r.h.s. of eq. (2.8) for X 2 (τ ) (see appendix C for details): where we have used eq. (3.8). The second piece in eq. (2.8) is easily computed as (see JHEP05(2016)008 Note that the terms linear in τ in the above small-time expansions cancel between the two contributions to X 2 (τ ) . After also subtracting the square X(τ ) 2 of the average value, we obtain the variance: Remarkably, the terms of O τ 2 have cancelled between X 2 (τ ) and X(τ ) 2 , so the net contribution to σ 2 ε (τ ) at small time is quartic in τ . This result is very interesting: it shows that, in the physically interesting regime at πτ 2 1, the dispersion σ ε (τ ) in the energy loss at large angles is comparable with the respective expectation value: σ ε (τ ) ∼ ε(τ ) ∼ τ 2 . A physical interpretation of this result will be presented in the next subsection. For latter purposes, we notice that, to order τ 4 , eq. (3.13) can be equivalently rewritten as ε 2 (τ ) (4/3) ε(τ ) 2 . For larger times πτ 2 1, on the other hand, the variance rapidly vanishes, as is should be expected: at large times, the whole energy is lost towards the 'condensate' (meaning, towards the medium), so the fluctuations vanish.
The results obtained in this section are illustrated in figure 5. As manifest there, the event-by-event fluctuations in the energy loss at large angles are numerically important (comparable to the respective average value, albeit somewhat smaller) for all times τ < 1, i.e. t < t br (E). Their relative strength, as measured by the ratio σ ε (τ )/ ε(τ ) , is particularly large at τ 0.5, which is the interesting regime for the phenomenology of highenergy jets at the LHC. For instance, for a realistic medium size L = 4 fm and a LP with E = 100 GeV, one finds (usingq = 1 GeV 2 /fm andᾱ s = 0.3) a value τ = L/t br (E) 0.3, for which the fluctuation to average ratio is σ ε / ε 0.5.

Physical discussion
An important result of the previous analysis is that, in the kinematical range of interest for di-jets at the LHC, the event-by-event fluctuations in the energy loss at large angles are comparable to the respective average value. Specifically (cf. eqs. (3.13) and (3.6), that we here rewrite in physical units), To understand this result, one needs to identify the physical origin of the various terms which occur in the small-τ expansion in eq. (3.11), notably the quartic term -the only one to survive in the final result for the variance, eq. (3.13).
To that purpose, one should recall that the double integral x , τ ) measures the total energy squared carried by pairs of gluons from the jet; hence, it is naturally biased towards large values for x and x .
By inspection of the calculations in appendix C, it is rather easy to trace back the term linear in τ , which is the dominant term at small τ . This is generated by a truly democratic branching of the LP: the latter splits into two daughter gluons with x 1 x 2 1/2, which are eventually measured: x x 1 and x x 2 (see figure 6 left). Such a hard democratic branching has a low probability, of order τ = t/t br (E) 1, but it generates a pair of gluons with the highest possible energies, hence it dominates the double integral in eq. (3.11). The order of magnitude of this contribution is set by the probability ∼ O(τ ) for the hard branching.
Furthermore, the term quadratic in τ comes from processes in which the LP radiates with probability of O(1) a primary gluon with energy fraction x 2 ∼ τ 2 (i.e. with energy ω 2 ∼ ω br (t)), which then initiates a democratic cascade (with probability of O(1), once again). In this case the measured gluons are the LP itself (with energy fraction x ∼ x 1 ∼ 1) and one of the relatively hard gluons in the democratic cascade generated by x 2 , namely a gluon having x ∼ x 2 ∼ τ 2 . Such a branching sequence is illustrated in figure 6 right. (The symmetric process, where x 1 and x 2 exchange their roles, gives an identical contribution.) The overall order of magnitude is now set by the total energy ∼ τ 2 carried by the softest measured gluon x .
Finally, the term of O τ 4 , which is the most interesting one for our purposes, receives contributions from the processes previously discussed in relation with eq. (3.10) and illus-JHEP05(2016)008 We see that the dominant contribution, of O τ 2 , to the dispersion σ ε (τ ) at small τ has the same physical origin as the respective contribution to the average energy loss ε(τ ) : the relatively hard primary emissions with x ∼ τ 2 , or ω ∼ ω br (t). This observation naturally explains the result in eq. (3.15), via the following chain of arguments: (a) the average number of such primary gluons is parametrically of order one, (b) their emissions are random and quasi-independent from each other, hence (c) the fluctuations in their number are of order one as well.
Point (a) follows from the very definition of the scale ω br (t), as the energy of the gluon emissions which occur with probability of order one during the time t. This is further corroborated by the fact that both the average energy loss, eq. (3.6), and the typical energy loss by the LP, cf. the discussion after eq. (3.1), can be expressed as ω br (t) times a number which is parametrically of O(1).
Point (b) can be understood as follows: the successive emissions of primary gluons with ω ∼ ω br (t) are uncorrelated with each other, since they cannot overlap in time. Indeed, these emissions are distributed in time during the interval t, the formation time for each of them t f (ω br (t)) ∼ᾱ s t is parametrically shorter than t, and the overall number of such emissions is of O(1). There is of course a correlation introduced by energy conservation, but this is rather weak, since the total energy carried by these emissions, of order ω br (t), is much smaller than the available energy, which is E.
Finally, point (c) is an immediate consequence of the two previous points. Together with point (a), it ultimately implies that both the average energy loss at large angles and its typical fluctuations should be of order ω br (t), in agreement with eq. (3.15).

JHEP05(2016)008
The above arguments also explain the gross structure of the LP peak in the gluon spectrum (3.1), i.e. the fact that the width δx ∼ τ 2 of this peak is of the same order as the typical energy loss by the LP, as measured by the shift 1 − x p ∼ τ 2 in the position of its maximum. The emphasis on the typical energy loss by the LP in the above discussion is indeed important, since this quantity is different from the respective average quantity ∆E : the former is associated with relatively soft primary emissions, with energies ω ω br (t), which occur event-by-event, whereas the latter is controlled by the hardest possible emissions, with energy ∼ ω c (t) ≡qt 2 [12,14]. Such hard emissions are rare events, which occur with a probability of O(ᾱ s ) (cf. eq. (2.2)), but they take away a large amount of energy and thus give the dominant contribution ∆E ∼ᾱ sq t 2 to the average energy lost by the LP.
The previous discussion also suggests that the emission of soft primary gluons should be a Poissonian process (see also ref. [26] for a similar argument). That is, the probability to radiate n gluons during a time t, such that each gluon is soft relative to its emitter (the LP), but harder than some 'infrared' scale ω 0 , is of the Poisson type, with the splitting rate given by eq. (2.2): This distribution predicts that the average number of primary emissions and its variance are equal to each other, which truly means that there are no correlations: n(n−1) = n 2 . Specifically, Notice the importance of introducing the 'infrared' scale ω 0 when studying the gluon distribution: without such a cutoff, the number of soft gluons would be infinite, as already obvious from eq. (2.2). By varying this cutoff ω 0 one can explore various regimes. In particular, eq. (3.17) confirms that both the average number of gluons with ω ∼ ω br (t) and its fluctuations are of order one. 9 It also predicts that the average number increase quite fast when decreasing ω 0 below ω br (t), while the relative strength of the fluctuations decreases: σ n / n = 1/ n 1 when ω 0 ω br (t). This last conclusion must be taken with a grain of salt, since the primary gluons with ω ω br (t) cannot be distinguished from the gluons with similar energies which are produced via multiple branching. The latter are expected to dominate the soft part of the gluon distribution and for them the typical fluctuations are parametrically large. This will be confirmed by the calculations in the next section.
Yet, there are quantities, including observables, which are sensitive to the primary gluons. One such an observable, as we have seen, is the energy loss at large angles; but this probes only to the relatively hard primary emissions, with ω ∼ ω br (t), so it cannot test the Poisson distribution (3.16) at softer energies. A better candidate in that sense is the energy JHEP05(2016)008 loss ∆E by the LP, which is by definition controlled by the primary emissions. As observed in [26], the precise distribution of ∆E (and not only its average value ∆E ∼ᾱ s ω c ) is important for understanding the quenching of the hadronic spectra at large p T . This distribution P(∆E) has been computed in ref. [26] under the assumption that primary emissions are Poissonian. Remarkably, the function P(∆E) thus obtained coincides with the shape of the LP peak in the gluon spectrum (3.1), with the identification ∆E ≡ E(1−x) and for 1 − x 1. We recall that the spectrum (3.1) has been obtained via a more general calculation, which includes multiple branchings and makes no special assumption about the primary gluons. In our opinion, this agreement between eq. (3.1) with x 1 and the result for P(∆E) in [26] provides a strong argument in favor of a Poisson distribution for the primary emissions. Further evidence in that sense will be presented in the next section.

Event-by-event fluctuations in the gluon distribution
At several places in the previous discussion, we made remarks about the correlations in the gluon distribution produced by the branching process. In section 2.2 we argued that gluons should be generally correlated with each other, because they have common ancestors. Later on, in section 3.3, we argued that the emissions of primary gluons should be independent of each other, because they have no overlap in time. In this section, we would like to substantiate such previous comments via explicit calculations of the gluon multiplicity event-by-event and in particular demonstrate that they are indeed consistent with each other: gluons can be mutually correlated, or uncorrelated, depending upon their production mechanism and upon their energies.
The random variable that we shall consider to that aim is the number N (t, x 0 ) of gluons with energy fractions larger than some infrared cutoff x 0 1 at time t, Via the correspondence between the low energy cutoff ω 0 ≡ x 0 E and the angular opening θ 0 of the jet, namely sin θ 0 Q s /ω 0 with Q 2 s (t) =qt, this quantity N (t, x 0 ) can also be viewed as the event-by-event number of gluons inside a jet with opening angle θ 0 . This looks similar to the jet energy fraction X(τ, x 0 ) introduced in eq. (3.7), but one should keep in mind that these two quantities probe very different aspects of the gluon distribution inside the jet (θ < θ 0 ): the energy fraction X(τ, x 0 ) is controlled by the hardest gluons and notably by the LP, whereas the gluon number N (t, x 0 ) is rather controlled by the softest gluons, with energies close to the infrared cutoff ω 0 . This will become clear after computing the first two moments of the gluon number distribution, that is, the mean number

JHEP05(2016)008
and its second moment, which also determines the variance (recall eq. (2.9)), The last equation can be rewritten as N 2 = N (N − 1) + N , where N (N − 1) is the second factorial moment (for gluons with x ≥ x 0 , of course). As before, these quantities will be computed for the simplified kernel K 0 , that is, by using eqs. (3.1) and (3.9) for the gluon spectrum and for the (energy-weighted) gluon pair density, respectively. As just mentioned, the gluon number distribution is very sensitive to the actual value of the infrared cutoff x 0 , hence the importance of properly chosing this value to match our physical purposes. Clearly, one needs x 0 1 in order to describe radiation. To probe multiple branching at small times τ 1, one needs the stronger condition x 0 τ 2 , that is, ω 0 ω br (t). For consistency with our approximation scheme, this cutoff ω 0 must be harder than the medium scale T , where the branching dynamics gets modified by elastic collisions and thermalization [42]; hence, we shall choose 10 T /E x 0 1. Note that we implicitly assume the separation of scales T ω br (t) E, which is indeed well satisfied in practice. We furthermore assume that the medium acts as a 'perfect sink' at the low energy end of a gluon cascade: it efficiently absorbs the energy carried by the soft branching products with ω T without modifying the branching dynamics at higher energies ω ≥ ω 0 . The recent study in ref. [42], which addressed the interplay between branching and thermalization, demonstrates that this assumption is rather well satisfied for sufficiently small times t t br (E), or τ 1.

The average gluon number
The calculation of the average gluon number according to eq. (4.2) is straightforward and gives where we have changed variables according to , and the approximations in the last line hold for πτ 2 1 and x 0 1. The above result for N (t, x 0 ) is illustrated in figure 7 for several values of τ ≤ 1.
The final result at small τ can be written as N (t, x 0 ) 1 + n(t, x 0 ) , where the 1 refers to the LP, while n(t, x 0 ) = 2τ / √ x 0 is the average number of radiated gluons with 10 Within the theory for turbulence, such an intermediate range in energies, which is well separated from both the 'source' and the 'sink', is known as the 'inertial range' [46,47]. x ≥ x 0 . The above integral yielding n(t, x 0 ) is dominated by its lower limit u = u 0 , or x = x 0 ; this shows that n(t, x 0 ) is essentially the same as the number of gluons with energy fraction x 0 . Hence, by studying the statistics of gluons with an infrared cutoff x 0 1, we actually describe the event-by-event distribution of the particles propagating at angles θ θ 0 , with θ 0 as defined after eq. (4.1).

JHEP05(2016)008
The above estimate for n(t, x 0 ) at small τ is identical with the respective prediction, eq. (3.17), of the Poisson distribution (3.16), which was supposed to refer to primary gluons only. That is, the average multiplicity of the radiated gluons with x ≥ x 0 is not modified by multiple branching -it is formally the same as produced via radiation by the LP alone. This is in line with the previous observation that the gluon spectrum at x 1, eq. (3.3), is a fixed point of the branching dynamics. Of course, the multiple branchings do lead to the copious production of soft gluons with x τ 2 , but these gluons decay as fast as they are produced and leave behind only very soft gluons, with energies ω T , which disappear in the plasma (formally, they accumulate into a condensate at x = 0). Such gluons are not counted in our estimate for N (t, x 0 ) , which is valid only for x 0 > T /E.
Another interesting consequence of the final result in eq. (4.4) is that the number n(t, ω 0 ) of soft radiated gluons is independent of the energy E of the LP (so long as ω 0 E and for sufficiently small times t t br (E)). This number becomes large, n(t, ω 0 ) > 1, when the cutoff is sufficiently soft, ω 0 ω br (t). We thus conclude that a high-energy jet with E ω br (L) produces a large number of soft gluons which propagate at large angles; this number is roughly independent of E and grows linearly with the distance L travelled by the jet through the medium.
Note finally that n(t, ω 0 ) is rapidly increasing when decreasing the energy cutoff ω 0 , or, equivalently, when increasing the jet opening angle θ 0 ∝ 1/ω 0 (see figure 7). This should be contrasted to the respective behavior of the jet energy, which is almost independent of JHEP05(2016)008 ω 0 (or θ 0 ) when ω 0 is small enough (recall figure 3). This difference is easy to understand: when decreasing ω 0 (or increasing θ 0 ), one includes new gluons within the jet, whose number is increasing as 1/ √ ω 0 (or equivalently like √ θ 0 ), but which carry only little energy ∝ √ ω 0 .

The gluon number fluctuations
We now turn to the gluon number correlations, as measured by the second factorial moment with D (2) (x, x , τ ) given by eq. (3.9). This calculation turns out to be more complicated and we have not succeeded in obtaining an exact result in closed form. But for the most interesting physical regimes, we have obtained accurate results with a transparent physical interpretation. In fact it is easy to anticipate these results via simple arguments that we now present. As before, we assume x 0 1 and focus on relatively small times τ 1, when the LP still exists.
The factorial moment N (N − 1) (τ, x 0 ) counts twice the number of gluon pairs in which both gluons have energy fractions larger than, or equal to, x 0 . Clearly, there are two types of such pairs: (A) those made with the LP plus a radiated gluon, and (B) those in which both gluons are radiation products. It is furthermore clear that the average number of pairs of types A is the same as the average number of radiated gluons with . Concerning the pairs of type B, we expect their number to be dominated by those pairs in which both gluons are soft, meaning that the integrals in eq. (4.5) are controlled by x, x 1. In this regime, one can use the simplified version of the pair density in eq. (3.10), to deduce where the integrals are indeed controlled by the lower limit (hence, the upper limit has been arbitrarily set to 1). Adding the previous results, one finds This turns out to be indeed the right result, at least in the two limiting situations where we have been able to approximately compute the integrals in eq. (4.5), namely (i) the multiple branching regime at x 0 τ 2 1 and (ii) the perturbative regime at τ 2 x 0 1 (see appendix D for details). The transition region around x 0 = τ 2 turns out to be more difficult to study, but the fact that we obtain identical results when approaching this region from both sides makes us confident that the transition is quite smooth and reasonably well described by eq. (4.7).
The r.h.s. of eq. (4.7) coincides with the respective prediction of the perturbative expansion in the number of branchings to second order: the first term, linear in τ , is the same as the result of a single emission by the LP, while the second term, quadratic in τ ,

JHEP05(2016)008
would be also obtained by summing the 3 processes which involve exactly 2 branchings, cf. figure 4. We thus conclude that the 2nd factorial moment is not modified by multiple branchings, similarly to the mean multiplicity. Once again, this feature is a non-trivial consequence of the wave turbulence, i.e. of the fact that the energy flux generated via democratic branchings is independent of x.
Using eq. (4.7) together with eq. (4.4), one can immediately compute the variance in the number of gluons with x ≥ x 0 : Note that, when rewritten in terms of the dimensionfull quantities t and ω 0 , both the second moment (4.7) and the variance (4.8) scale like functions of ω br (t)/ω 0 ∝ t 2 /ω 0 and are independent of the energy E of the LP. A similar behavior has been observed for the average multiplicity (4.4).
The physical consequences of these results will be now separately discussed for the two limiting regimes aforementioned.
Besides confirming the previous estimate in eq. (4.7), this result also specifies the respective error. 11 Clearly, in this regime, the term quadratic in τ in the r.h.s. of eq. (4.9) dominates over the linear one. By keeping only this dominant term, together with the corresponding approximation N 2τ / √ x 0 1 for the average gluon number, one can write This shows that the second factorial cumulant N (N − 1) − N 2 , which we recall is a measure of the correlations, is parametrically as large as the 'disconnected' 2-point function N 2 (and much larger than the mean number N ). This means that correlations are large.
Recalling the physical origin of the term quadratic in τ , cf. eq. (4.6), it is clear that the net correlation, namely The calculations in appendix D also confirm the physical origin of the 2 terms in the r.h.s. of eq. (4.9): the quadratic term is controlled by very soft gluons, with x ∼ x ∼ x0; the linear term arises by integrating over x ∼ x0 and 1 − x πτ 2 , or vice-versa; that is, this term measures the pairs made with one soft gluon and the LP.

JHEP05(2016)008
comes from processes in which the two measured gluons belong to a same mini-jet, i.e. they have a soft common ancestor, cf. figure 1 left. The contributions to N (N − 1) coming from two different mini-jets, cf. figure 1 right, cancel against the disconnected piece N 2 . Soft gluons which belong to different mini-jets are uncorrelated with each other, because so are the primary gluons which have initiated those mini-jets in the first place (see section 4.2.2 below).
To the accuracy of interest, the variance can be read from either eq. (4.11), or eq. (4.8); one thus finds σ 2 N (τ, x 0 ) 2τ 2 /x 0 . This is large, σ 2 N N 1, showing that the gluon distribution at small x is overdispersed (as compared to the Poisson distribution, which would predict σ 2 N = N ). The above results can be summarized as when t t br (E) and ω 0 ω br (t) , (4.12) where we have restored the physical units to emphasize that, to the accuracy of interest, the result is independent of E. Eq. (4.12) is similar to eq. (3. 15), in that it shows that the typical fluctuations in the gluon number distribution at small x are as large as the respective average value. We expect factorization properties similar to eq. (4.10) also for the higher moments: where the factors C p (known as 'reduced moments') are pure numbers. Indeed, as shown by the example of the 2-point function (3.10), it is natural to expect the factorization of the p-point correlation function D (p) (x 1 , x 2 , . . . , x p ) for small values of its arguments, x i 1: the soft gluons evolve independently from each other after the splitting of their last common ancestor. Notwithstanding, the correlations are strong, because the information about the LCA cannot be lost: it determines the topology of the cascade.
A probability distribution for which the (factorial) moments of order p scale like the p-th power of the mean number is said to obey KNO scaling [43]. This is tantamount to saying that the product N P N between the mean number and the respective probability distribution P N depends upon the random variable N and the various parameters of the distribution (in our case τ and x 0 ) only via the ratio N/ N . This has interesting consequences for the phenomenology: the ratio N 2 / N 2 is predicted to be independent of τ and x 0 (and similarly for the higher reduced moments C p ). This could be explicitly checked in the data -say, in the asymmetric di-jet events in AA collisions -, by measuring the event-by-event distribution of soft gluons propagating at large angles.
One additional reason why the KNO scaling looks appealing in the present context is because a similar scaling (albeit with a different value for the second reduced moment C 2 ; see below) has been observed to hold for jets evolving in the vacuum [44]. In that case, the branchings are triggered by the virtuality of the leading particle, so the branching probability is different -it is given by the DGLAP splitting functions. This difference has indeed consequences for the details of the distribution, like the values of the reduced JHEP05(2016)008 moments C p . For jets evolving in the vacuum, one finds C 2 = 4/3 in the leading doublelogarithmic approximation [44]. That is, the vacuum analog of our eq. (4.10) reads 12 N (N − 1) = 4 3 N 2 (jets in the vacuum). (4.14) This vacuum distribution looks already quite broad (say, as compared to a Poisson distribution, for which N (N −1) = N 2 ), but the one generated via medium-induced branchings, for which C 2 = 3/2, is even broader. So, the fluctuations in the gluon distribution produced by a jet propagating through a dense medium should look even larger (in appropriate units, cf. eq. (4.13)) than for a jet which propagates in the vacuum.

τ 2
x 0 1: primary gluons and the Poisson distribution In this regime, gluon branchings have a low probability, hence the average number of radiated gluons is small, n(t, x 0 ) 2τ / √ x 0 1, and fully controlled by direct emissions by the LP: these are primary gluons. Concerning the second moment, the calculations in appendix D yield (cf. eq. (D.23)) where the linear (quadratic) term represents the leading (subleading) contribution. The perturbative expansion in the number of branchings is now meaningful and the two terms in the r.h.s. of eq. (4.15) can be physically (and not only formally) recognized as the results of a single emission (for the term linear in τ ) and, respectively, of a double emission (for the quadratic term). The variance too is dominated by the linear term in the r.h.s. of eq. (4.8) and hence it is equal to the average number n(t, x 0 ) 2τ / √ x 0 of primary gluons. This is consistent with our conclusion in section 3.3 that the primary emissions with x 1 should obey a Poisson distribution, cf. eqs. (3.16)- (3.17). This argument can be made more precise, as we now explain.
To that aim, notice that the quadratic term in the r.h.s. of eq. (4.15) receives 2/3 of his strength from processes involving the emission of 2 primary gluons, cf. figure 4. Hence, after subtracting 1/3 of this term (which is associated with the production of 2 secondary gluons, as shown by the first process in figure 4), the remaining contribution, that is, 12 An elementary probability distribution which is known to lead to KNO scaling in the regime of large numbers is the negative binomial distribution (NBD), whose definition involves a parameter r. The relations (4.10) and (4.14) are consistent with the NBDs with r = 2 and r = 3, respectively. Yet, we see no physical reason why the NBD with r = 2 should truly apply to the gluon distribution produced via medium-induced branchings. Indeed, for a jet decaying in the vacuum, for which all the factorial moments can be explicitly computed, the NBD with r = 3 fails to describe these moments for large values of pthat is, it does not reproduces the respective values of Cp for p ≥ 3 [44].

JHEP05(2016)008
involves only primary emissions. For such emissions, it makes sense to decompose the random variable N as N = 1 + n (the LP plus the number n of primary gluons). We then easily deduce 17) or, equivalently, σ 2 n = n , which is the hallmark of the Poisson distribution, cf. eq. (3.17).

Summary and physical discussion
In this paper, we have for the first time investigated the effects of event-by-event fluctuations in the gluon distribution produced via medium-induced branching by an energetic jet propagating through a dense QCD plasma. We have identified a characteristic pattern for the medium-induced radiation, which can be observed on an event-by-event basis, albeit the fluctuations from one event to another are predicted to be large. This pattern is characterized by the production of many soft gluons -most of them, as soft as the medium scale T -which propagate at large angles w.r.t. the jet axis and collectively carry a large amount of energy. Importantly, this soft medium-induced radiation is predicted to be geometrically separated -in the sense of propagating at larger angles -from the vacuum-like radiation associated with the virtualities of the jet constituents. Accordingly, the physical effects that we found, like the importance of fluctuations, should be also visible in the actual experiments, where both sources of radiation coexist with each other: the medium-induced radiation is expected to control the (energy and hadron number) distribution produced by the jet at large angles. This being said, the vacuum-like fragmentation is an essential ingredient of the overall picture and must be included in realistic applications to phenomenology. This is itself a random process, hence it should amplify the fluctuations and also the inmedium energy loss. Indeed, any of the (relatively hard) partons produced via vacuum-like fragmentation that can be resolved by the medium will play the same role as the 'leading particle' in our discussion -it will act as a source for soft medium-induced radiation, thus enhancing the energy loss at large angles (see e.g. the discussions in [33,45]).
The overall picture depends upon the ratio between the energy E of the leading particle and the medium scale ω br (L) =ᾱ 2 sq L 2 (the characteristic energy for multiple branching), or, equivalently, between the branching time t br (E) = (1/ᾱ s ) E/q for the LP and the distance L travelled by the jet through the medium. For not so high energies, such that E ω br (L) (or t br (E) L), the LP disappears via democratic branching and its whole energy is ultimately transmitted to soft gluons propagating at large angles. In the highenergy regime at E ω br (L) (or t br (E) L) -the most interesting regime for the phenomenology of di-jet asymmetry at the LHC -the LP survives in the final state and the energy transferred at large angles in a typical event is independent of E and of the order of the medium scale ω br (L). Moreover, the event-by-event fluctuations in the energy loss at large angles are of order ω br (L) as well -that is, they are as large as the respective average value.

JHEP05(2016)008
The physical interpretation of this result is quite interesting: the energy gets transmitted to softer and softer quanta via successive quasi-democratic branchings. The energy flux generated by such branchings is quasi independent of ω. Hence the energy loss at large angles is controlled by the hardest partons which can initiate democratic cascades during a time ∼ L -those with energies ω ∼ ω br (L). Such gluons are emitted by the LP with a probability of order one and their emissions are quasi-independent of each other. Hence, both the average number of such gluons and its typical fluctuations are of O(1). In turn, this implies that both the average energy loss at large angles and the respective dispersion should be of order ω br (t), which is what we found (cf. eq. (3.15)).
Another remarkable feature of the gluon distribution x(dN/dx) produced by an energetic jet with E ω br (L) is the fact that, for the soft gluons with x = ω/E 1, this depends only upon the dimensionless ratio τ / √ x = [ω br (L)/ω] 1/2 . In particular, it is independent of E. This property too is a consequence of wave turbulence, i.e. of the fact that the gluon distribution at small x is not modified by multiple branching. For instance, the gluon spectrum D(x) with x 1 has the scaling form visible in eq. (3.3), while the 2-point function D (2) (x, x ) with x, x 1 factorizes as shown in eq. (3.10). A similar factorization property is expected for the higher-point correlations. It reflects the fact that soft gluons evolve independently from each other after the splitting of their last common ancestor.
These special properties of the gluon distribution at small x have interesting consequences for the event-by-event number N (ω 0 ) of gluons with energies ω ≥ ω 0 . First, the multiplicity distribution too is independent of E and solely a function of ω br (L)/ω 0 . Second, in the non-perturbative regime at ω 0 ω br (L), where the multiplicity is large, the higher moments N p (ω 0 ) are predicted to obey KNO scaling. For instance, the ratio C 2 ≡ N 2 / N 2 should be a pure number, independent of all the physical parameters of the problem (like E, L,q, or the infrared cutoff ω 0 ). This strong prediction should be taken with a grain of salt, as it assumes that the only source of fluctuations is the stochasticity of the medium-induced branching process. In the actual collisions though, there are also other stochastic aspects, like the vacuum-like fragmentation and the geometry of the hard process giving rise to the leading particle. Hence, the distance L travelled by the jet trough the medium is itself a random variable and equations like (4.10) or (3.13) should be more properly read as where the various averages also include the event-by-event distribution of L. Such relations are both interesting and non-trivial: not only they are still independent upon physical parameters like E,q, or ω 0 , but they also show that, by combining event-by-event measurements of the soft hadron distribution at large angles with the results of the present theoretical analysis, one can better constrain the distribution of L in the experiments. The existence of large event-by-event fluctuations in the energy loss at large angles should also have important consequences for the phenomenology of di-jet asymmetry, as recently demonstrated in ref. [45], via Monte-Carlo simulations using the event generator JEWEL [48]. Namely, the fluctuations can lead to large values for the di-jet asymmetry A J event-by-event, including in events where the in-medium path lengths of the 2 jets are JHEP05(2016)008 quite similar. In other terms, due to the large fluctuations in the energy loss, one loses the generally expected relation between the distribution of A J and that of the difference L 2 − L 1 between the distances, L 1 and L 2 , travelled by the 2 jets through the medium.
To explain this, let us first recall that the di-jet asymmetry is expressed in terms of the energy imbalance 13 , (5.2) between the leading jet (E 1 ) and the subleading one (E 2 ), with E 1 > E 2 due to the 'trigger bias' (i.e. the fact one enforces different kinematical constraints on the two jets). In the equation above, One generally assumes that the difference E 2 − E 1 in the energy losses is controlled by the difference L 2 − L 1 between the path lengths, but in reality this correlation can be loosen by fluctuations. For instance, if one assumes a class of events with fixed values for L 1 and L 2 (this is experimentally unrealistic, but can be artificially simulated via Monte-Carlo generators), the value of E 2 −E 1 can still differ from event to event, due to fluctuations in the (vacuum-like or medium-induced) radiation by each of the two jets. The main numerical observation in [45] is that, in central collisions at least, the path-length difference L 2 − L 1 is typically small, yet the event-by-event asymmetry A J is generally large (and only weakly correlated with L 2 − L 1 ), due to large fluctuations in the radiative processes. Our present analysis corroborates and strengthens the conclusions in [45] by clarifying the microscopic origin of the fluctuations in the medium-induced radiation and the associated energy loss.
To more explicitly demonstrate the role of fluctuations for the di-jet asymmetry, let us consider the fluctuations in the energy loss, as computed in section 3. For convenience, we replace A J by the simpler quantity 14 ε 2 − ε 1 and give up the kinematical cuts which distinguish between the leading and the subleading jet -so that the difference ε 2 − ε 1 can have either sign. Its average value ε 2 − ε 1 ∝ L 2 2 − L 2 1 (cf. eq. (3.5)) is a measure of the difference between the path-lengths and thus it would vanish in a Gedanken (Monte-Carlo) 'experiment' where one fixes L 1 = L 2 event-by-event. However, without the kinematical cuts, the di-jet asymmetry is better characterized by the respective dispersion, easily computed as where the last estimate follows from eq. (3.13). Notice that the r.h.s. of eq. (5.3) is proportional to the sum, L 2 1 + L 2 2 , of the average squared path-lengths. Hence, even if one 13 In the experimental definition of AJ , the total jet energy Ei is replaced by its transverse momentum p ⊥i , i.e. its total momentum in the plane perpendicular to the collision axis (see e.g. [1,2]). For jets propagating at nearly central rapidities, one has p ⊥i Ei. For the sake of illustrating our argument here, it is more convenient to use the jet energies, instead of their transverse momenta.

JHEP05(2016)008
artificially fixes L 1 = L 2 , one still finds a non-zero di-jet asymmetry A J ∝ (ε 2 − ε 1 ) 2 , which is moreover large -of the order of the average fractional energy loss ε by any of the two jets -and which in particular scales with the path-length like L 2 .
Of course, in actual experiments, the path-lengths L 1 and L 2 are randomly distributed, but the above argument suggests that the contribution to A J coming from fluctuations in the medium-induced radiation can be as large as the respective contribution from differences in path-lengths -perhaps even larger in the case of central collisions, where the path lengths for the 2 jets are not very different from each other [45].
We conclude this physical discussion with a few numerical estimates, to give a flavor about the size of the fluctuations in physical units and for phenomenologically interesting situations. We consider a leading particle with E = 100 GeV which propagates through the medium along a distance L for which we shall consider 3 representative values: L = {2; 4; 6} fm. As before, we choseᾱ s = 0.3 andq = 1 GeV 2 /fm. With these choices, one finds t br (E) 15 fm, hence the 'high-energy' condition t br (E) L is reasonably well satisfied for all the interesting values of L. Corresponding to the 3 values for L aforementioned, we find: ω br (L) {2; 8; 18} GeV. Using eq. (3.15), the average energy loss at large angles and the associated dispersion are roughly estimated as E(L) {6; 25; 57} GeV and, respectively, σ E (L) {3; 15; 33} GeV. This shows that, e.g. for L = 4 fm, the expected range for the variation of the event-by-event energy loss is 10 < E < 40 GeV.
Consider similarly the statistics of the number of gluons with energy ω ≥ ω 0 , with ω 0 chosen of the order of the medium temperature T , say ω 0 = 0.5 GeV. The mean number of radiated gluons for the 3 choices of the medium size L is then estimated as n(L, ω 0 ) = 2 [ω br (L)/ω 0 ] 1/2 {4; 8; 12}. Note that this number grows linearly with L. As explained after eq. (4.4), most of these gluons are soft, say with energies below 2 GeV, hence they propagate at large angles, outside the reconstructed jet in the final state. The variance is estimated according to eq. (4.8) as σ 2 N (L, ω 0 ) {12; 40; 84}. Hence, for L = 4 fm we expect n(ω 0 ) to take values within the range 2 < n(ω 0 ) < 14. In an actual experiment, these numbers should further increase due to the gluon fragmentation into hadrons.
Quite remarkably, the above numbers agree quite well with the phenomenology of dijet asymmetry at the LHC [1][2][3][4][5][6][7][8][9]. Such an agreement should be taken with a grain of salt, since our present description is rather idealized: we have ignored several important dynamical ingredients, like the elastic collisions off the medium constituents, the longitudinal expansion of the medium, the virtualities of the partons in the jet, the gluon fragmentation into hadrons, and the additional sources of fluctuations, besides the randomness of the branching process. It is in principle possible to systematically include such effects within Monte-Carlo event generators, such as MARTINI [29] or JEWEL [48], and it is of outmost importance to do so in practice, in order to promote the physical picture developed in this paper into a realistic approach to the phenomenology.

JHEP05(2016)008
and Yacine Mehtar-Tani during the early stages of this work. This work is supported by the European Research Council under the Advanced Investigator Grant ERC-AD-267258 and by the Agence Nationale de la Recherche project # 11-BS04-015-01.

A Master equations for the Markovian branching process
In this appendix, we shall construct the transport equations (2.13) and (2.14) obeyed by the energy density D(x, t) and by the gluon pair density D (2) (x, x , t), respectively. To that aim, we start with the master equations obeyed by the probability densities P n (x 1 , x 2 , · · · , x n |t) introduced in section 2.2. These equations can be easily established via elementary probabilistic considerations. We recall that P n (x 1 , x 2 , · · · , x n |t) is the semi-inclusive probability density for having a state with n gluons with energy fractions x i ≥ (i = 1, . . . , n), together with an arbitrary number of softer gluons with x < , which are not measured. The dependence upon the infrared cutoff will be generally omitted at intermediate steps since the final equations for the correlations admit a well-defined limit → 0.
In one evolution step, t → t + dt, this probability P n can increase due to the splitting of a gluon from a preexisting state with n − 1 gluons, but it can also decrease via the decay of one of the n gluons explicitly included in P n . The splitting rate can be read off eq. (2.12), conveniently rewritten as where we recall that x = ω/E is the energy fraction of the parent gluon and z is the splitting fraction; that is, the 2 daughter gluons have energies zω and (1 − z)ω, respectively. These considerations motivate the following evolution law for P n , involving a loss term and a gain term: dz K(z, x i ) P n (x 1 , x 2 , · · · , x n |τ ) Recalling equations (2.12) and (A.1), it is quite clear that the integral over z in the negative, 'loss', term has endpoint singularities at z = 0 and z = 1. These are regulated by the conditions that both daughter gluons be harder than the infrared cutoff: zx > and (1 − z)x > . Given the master equation (A.2) and the rule (2.10) for computing the expectation value O(τ ) of a generic observable, it is straightforward (albeit a bit tedious) to deduce an evolution equation for the latter. This is conveniently written as

JHEP05(2016)008
with the notations {x} = (x 1 , x 2 , . . . x n ) and The subscript i in eq. (A.4) refers to the gluon x i which splits in the evolution step under consideration. Note that, in general, (A.3) is not a closed equation, that is, the r.h.s. cannot be written as a linear operator acting on O(τ ) . As a first example, consider the gluon spectrum D(x, τ ), as defined in eq. (2.9). In this case, D n ({x}; x) = x n i δ(x i − x) and Note that, in this case, the r.h.s. of the above equation is "local in i" -that is, it depends upon the n-gluon configuration {x} = (x 1 , x 2 , . . . x n ) only via the energy fraction x i of the gluon which has split in this particular evolution step. Accordingly, after inserting eq. (A.5) into the general equation (A.3), one finds a closed equation for the gluon spectrum, namely In the first 'gain' term, the integral over z is restricted to z ≥ x [by the support of D(x/z, τ )] and has an endpoint singularity at z = 1. In the second 'gain' term, we similarly have 1 − z ≥ x and a singularity at z = 0. Finally, the negative, 'loss', term has two endpoint singularities, at z = 0 and z = 1. Each of these terms could be individually regulated by the infrared cutoff , but this is actually not needed: the singularities mutually cancel in the sum of the three terms, so the complete equation admits a well-defined limit → 0, as anticipated. Using the symmetry property K(z, x) = K(1 − z, x), it is clear that the two 'gain' terms give identical contributions. As for the 'loss' terms, this can be rewritten in a more convenient form by using K(z, x) = K(z)/2 √ x together with the identity Putting things together, we recognize eq. (2.13), as expected. The evolution of the gluon pair density D (2) (x, x , τ ) could be similarly treated, but the corresponding argument is considerably more tedious. An alternative formulation, which is more convenient in practice, relies on the following generating functional with u(x) an arbitrary 'source' function with support at 0 ≤ x ≤ 1. Probability conservation requires Z τ [u = 1] = 1. It is easily checked that the correlation functions of the gluon density can be obtained as functional derivatives of Z τ [u] evaluated at u(x) = 1. For instance Clearly, the generating functional can be viewed as a special kind of 'observable', hence it obeys an evolution equation with the structure shown in eq. (A.3). Using the latter together with simple manipulations, one finds (see also [23] for a more general equation of this type, which includes the effects of transverse diffusion) This is a linear, and closed, functional differential equation. The functional derivative δZ τ [u]/δu(x) plays the role of an 'annihilation operator' (it reduces by one the number of factors of u), the term quadratic in u describes the 'gain' effect associated with the branching process x → (zx, (1 − z)x), whereas the negative term linear in u is the corresponding 'loss' effect. By taking one functional derivative w.r.t. u(x) in eq. (A.10) and evaluating the result at u(x) = 1, one immediately recovers eq. (A.6) for D(x, τ ). By taking two such derivatives, one finds an equation for D (2) (x, x , τ ) which after simple manipulations takes the form shown in eq. (2.14).

B Computing the gluon pair density
In this appendix we shall derive eq. (3.9) for the gluon pair density. Our starting point is eq. (2.15), that we shall compute for the simplified kernel K 0 (z) ≡ 1/[z(1 − z)] 3/2 . Accordingly, the Green's function D(x, τ ) is given by eq. (3.1), while the source term takes the form (recall eq. (2.14)) (B.1) Using this together with simple algebraic manipulations, eq. (2.15) can be rewritten as We focus first on the integration in x 2 and denote After the change of variables u = x 2 −x 1−x 1 −x 2 , we see that the only combination of energy fractions that the integral is sensitive to is (1 − x 1 − x ) : The integration in u can be done analytically, by using

JHEP05(2016)008
Differentiating the above with respect to a, we obtain another useful relation, By combining the two above equations, we deduce and therefore We are left with the integration in x 1 , which now reads (B.9) This integral over x 1 is very similar to the one over x 2 that we have just performed, cf. eq. (B.3): The last integral in τ is easy to do by using dx x e −x 2 = − e −x 2 2 and yields the result in eq. (3.9).

C Computing energy fluctuations
In this appendix we present in some detail the calculations leading to eqs. (3.11) and (3.12). For convenience, we introduce the notations and we start by computing X 2 (τ ) A using eq. (3.9) for D (2) (x, x , τ ) : The integration in x can be simplified by making the change of variable u = x 1−x−x : where

JHEP05(2016)008
Using the result in eq. (B.5) we see that dI dx = − π x e −x which implies and therefore The remaining integral in x becomes simpler if one first takes a derivative w.r.t. τ : This integral is very similar to the one we have performed over x and gives By also using dx erf(x) = x erf(x) + e −x 2 √ π we get (C.9) which is equivalent to eq. (3.11). Now we move to the computation of X 2 (τ ) B . Using eq. (3.1) for D(x, τ ), this is rewritten as The last integral can be done using integration by parts and the result is (C.14) Combining this with eq. (C.11) gives the result in eq. (3.12).

JHEP05(2016)008
In order to get more physical intuition about the physical meaning of these results, it is interesting to see what are the integration regions which contribute up to O τ 4 in the small time regime at τ 1, or πτ 2 1. Our purpose is to substantiate the physical discussion in section 3.3.
It is easy to check that the leading and subleading contributions to X 2 (τ ) A , of O(τ ) and respectively O τ 2 , come from integrating over x and x with 1 − x − x ∼ πτ 2 . Indeed, values such that 1 − x − x πτ 2 are strongly suppressed by the exponentials in eq. (C.2). In the region where 1 − x − x πτ 2 , we can expand the exponentials in eq. (C.2) and then we are led to integrations of the type with n ≥ 1. (Notice that the term with n = 0 has cancelled between the two exponentials.) But clearly, these integrals have an endpoint singularity at 1 − x − x = 0 (for any n ≥ 1), which in the original integral in eq. (C.2) was cut off by the exponentials at values x and x such that 1 − x − x ∼ πτ 2 . Then the dominant contribution comes precisely from the would-be singular endpoint, that is, from the pairs (x, x ) with 1 − x − x ∼ πτ 2 .
To more finely disentangle contributions of various orders in τ , let us now have a look at the integral representation in eq. (C.6). The leading-order contribution, linear in τ , can be obtained by expanding the error functions inside the integrand to linear order in their arguments; this yields X 2 (τ ) A πτ 2 and also shows that the relevant values of x obey x ∼ 1/2. Together with 1 − x − x ∼ πτ 2 , this implies that the term of O(τ ) is generated by x ∼ x ∼ 1/2, as anticipated in section 3.3. Beyond linear order, the integrand in eq. (C.6) cannot be expanded anymore, as this would generate non-integrable singularities at x = 1. This shows that the contribution of O τ 2 comes from 1 − x ∼ πτ 2 , which together with 1 − x − x ∼ πτ 2 implies that this quadratic piece comes from pairs (x, x ) made with the LP plus one soft particle: x 1 and x ∼ τ 2 , or vice-versa. It is furthermore easy to check that the respective contributions to X 2 (τ ) B have a similar origin: the contribution to eq. (C.10) of O(τ ) comes from x ∼ 1/2, whereas that of O τ 2 comes from x ∼ τ 2 .
It seems difficult to precisely characterize the integration region in x and x which contributes to X 2 (τ ) A at O τ 4 : this is a subsubleading piece, hence already the regions previously identified can give (subleading) contributions of O τ 4 . It is however clear that, to this order, one opens a new physical channel, which refers to the production of pairs (x, x ) where both particles are relatively soft, so that they are produced with a probability of order one, yet at the same time they are hard enough to contribute to the energyweighted pair density: these conditions select x ∼ x ∼ τ 2 . To see this, let us focus on the contribution of the soft modes to eq. (C.2). This can be estimated by introducing an upper cutoff x Λ ∼ τ 2 in the integrals over x and x , while at the same time keeping only the first non-trivial term, of O τ 2 , in the expansion of the exponentials (this expansion makes sense in the presence of the cutoff, since 1 − x − x 1 τ 2 ). One thus finds, to JHEP05(2016)008 parametric accuracy, 16) which is indeed of O τ 4 for the physically motivated choice x Λ ∼ τ 2 .

D Computing multiplicity fluctuations
In this appendix, we shall compute the second factorial moment of the gluon number distribution, N (N − 1) (τ, x 0 ), by explicitly performing the integrations over x and x which appear in eq. (4.5). After also using eq. (3.9) for D (2) (x, x , τ ), the starting point integral reads The integral over x can be simplified via the change of variables u = x 1−x−x , which gives where u 0 = x 0 1−x−x 0 . The integral over u is similar to that previously computed in eq. (C.13), hence After substituting the explicit form for the function G(a) from eq. (C.14), this result can be decomposed into two pieces, corresponding to the 2 terms in the r.h.s. of eq. (C.14). Specifically, N (N − 1) = N (N − 1) A + N (N − 1) B , where 1−x−x 0 (D. 4) and respectively In order to understand better the physical origin of the different terms it is useful to note that where the new integration variable t is related to the previous variables u and x appearing in (D.2) and respectively (D.1) (say, for the first exponential within the integrand there), via (D.7)

JHEP05(2016)008
(For the second exponential within the integrand one should simply replace τ → 2τ .) The first integral in eq. (D.6) corresponds to the first term in eq. (C.14) and is dominated by t ∼ √ a; this implies that the integral over x that led to the result in eq. (D.4) is controlled by its lower limit x ∼ x 0 . Similarly, the second integral in eq. (D.6) corresponds to the second term in eq. (C.14) and is controlled by t ∼ 1; accordingly, the integral over x leading to eq. (D.5) is controlled by x 1−x−x ∼ 1−x πτ 2 . We have not been able to analytically compute the remaining integrals over x, but we managed to obtain good approximations in two interesting limits, to be described in what follows.
D. 1 The case x 0 πτ 2 1 We start by computing N (N −1) A according to eq. (D.4). The region at 1−x−x 0 πτ 2 being exponentially suppressed, it is enough to consider the domain 1 − x − x 0 πτ 2 , meaning 1 − x πτ 2 + x 0 πτ 2 , where we have used the fact that πτ 2 x 0 . In other terms, we can substitute 1 − x − x 0 1 − x up to higher order terms, so the integration becomes After the change of variables u = x 1−x , the integral is similar to the one in eq. (C.13), so we have The term quadratic in τ / √ x 0 , which is the dominant term in this regime, comes from the lower limit x ∼ x 0 . The subleading term, linear in τ / √ x 0 , arises by integrating over x 1−x ∼ 1 πτ 2 , which for τ 1 is tantamount to x ∼ 1 − πτ 2 . Recalling that the whole contribution to N (N − 1) A comes from x ∼ x 0 , we conclude that the dominant term 6τ 2 /x 0 is generated by the region where both measured gluons live close to the infrared cutoff, x ∼ x ∼ x 0 , whereas the subleading term 2τ / √ x 0 comes from the pairs made with the LP plus a soft gluon: x 1 − πτ 2 and x ∼ x 0 , or vice-versa. Consider similarly the other contribution, N (N − 1) B in eq. (D.5). The terms involving error functions give only subleading contributions: indeed, the arguments of the error functions are much smaller than one except in the region where the integrand is already exponentially suppressed. The pole of the integrand at x = 1 is cut off by the exponential at 1 − x ∼ πτ 2 , hence its contribution to the final result is subleading too. We conclude that the dominant contribution comes from the lower limit x = x 0 , due to the 1/x 3/2 pole. So, we can write