Gluon density fluctuations in dilute hadrons

Motivated by the relation existing between the gluon density in a hadron and the multiplicity of the particles measured in the final state of hadron-nucleus collisions, we study systematically the fluctuations of the gluon density in onia, which are the simplest dilute hadrons, of different sizes and at various rapidities. We argue that the small and the large-multiplicity tails of the gluon distributions present universal features, which should translate into properties of the multiplicity of the particles measured in the final state of high-energy proton-nucleus collisions, or of deep-inelastic scattering at a future electron-ion collider. We propose simple physical pictures of the rare events populating the tails of the multiplicity distribution that allow us to derive analytical formulas describing these universal behaviors, and we compare them to the results of Monte Carlo simulations.


Introduction
Measurements of the number of particles produced in the final states of high-energy hadronnucleus (p-A) collisions have opened new windows in the study of dense partonic systems. Indeed, the multiplicity of particles detected in p-A collisions in the region of fragmentation of the hadrons is expected to be sensitive to the properties of their partonic content at the time of interaction. Nowadays, distributions of particle multiplicities are measured with great accuracy at particle colliders, most notably, in p-Pb collisions at the Large Hadron Collider (LHC) [1,2,3], and the possibility of relating the experimental data to the small-x dynamics of the parton distributions has triggered a lot of activity in the theoretical community [4,5,6]. A striking observation made in experiment is that multiplicity distributions in p-Pb collisions present a high-multiplicity tail that is much longer than in Pb-Pb collisions, as expected from the fact that p-A collisions are more fluctuation-dominated [7]. This tail may, then, contain valuable information about the wave functions of the projectile protons, and their event-byevent fluctuations.
Models of multiplicity fluctuations in the framework of high-energy hadronic collisions are mainly of two kinds. On the one hand, we have purely phenomenological models, that serve as initial conditions for hydrodynamic calculations, and that are typically based on rather ad hoc modifications of the Glauber model [8,9]. These models have largely grown in complexity over the past few years, and include now prescriptions to model the sub-nucleonic degrees of freedom in the proton [10,11,12,13]. Although such prescriptions are very successful in reproducing the experimental data, providing insight about the underlying dynamics of the parton distributions is beyond their scope. A different kind of calculations of particle multiplicities, that takes as input the configurations of gluons inside protons and nuclei, have instead been achieved within the color glass condensate picture [14,15,16,17,18] of high-energy quantum chromodynamics (QCD). These calculations are also very successful in phenomenological applications. A famous example is that of the number of gluons produced from the decay of color flux tubes in the glasma framework [19], which is distributed according to a negative binomial distribution, whose long tail at large multiplicity provides naturally a good description of proton-proton data.
In this paper, we work in the theoretical framework of high-energy QCD, and, following the picture of particle production introduced in Refs. [20,21], we argue that the multiplicity of particles probed around some off-forward rapidity in the region of fragmentation of the protons reflects, in each event, the integrated gluon density in the corresponding realization of the Fock state of the hadron at the time of interaction. In this picture, then, fluctuations of the multiplicity of particles are strictly related to the event-by-event fluctuations of the gluon density. Instead of studying this phenomenon directly in the case of proton-nucleus collisions at the LHC, which will require to introduce some amount of modelization for the evolution of the proton, we focus on the simpler case of onium-nucleus collisions, where one can gain a very solid theoretical understanding in controlled asymptotic limits that allow to study multiplicity fluctuations analytically. To this purpose, we shall work within the color dipole model (supplemented with an infrared cutoff for parton confinement), whose formulation is particularly well-suited to address this problem.
Our starting point is the main result of Ref. [21], namely, an analytical estimate of the behavior of the high-multiplicity tail of the gluon number density in a boosted onium. After reviewing this calculation, we extend the analysis to the opposite case of low-multiplicity events, and we derive a new formula for the behavior of the low-multiplicity region of the gluon number density. However, both our calculation and that of Ref. [21] rely largely on conjectures, and leave important parameters undetermined. The main thrust of this paper is, eventually, that of checking the validity of these derivations, and better understand them, by means of extensive Monte Carlo simulations of the small-x evolution of the Fock state of an onium in the dipole picture. Our goal is that of showing that the asymptotics of gluon density fluctuations in an onium present robust universal features.
Our motivation for pursuing theoretical studies of multiplicity fluctuations in the dipole model, pioneered over 20 years ago by Salam [22], is twofold. First, the theoretical understanding of dipole evolution has improved since then, as well as the numerical capabilities, making possible much more accurate evaluations of distributions. Second, and more importantly, there is now a strong motivation to better understand this physics, since the LHC is taking data for which such studies are relevant and timely. Our work may also be of interest for a future electron-ion collider, and actually, easier to connect to the experimental data in that case: Indeed, in a high-energy scattering, the interaction between an electron and an ion is mediated by a photon, whose qq (onium) component of the wave function is perfectly determined in the framework of quantum electrodynamics (QED).
Our paper is organized as follows. In Sec. 2, we recall the connection between the final-state multiplicity in hadron-nucleus collisions and the integrated gluon number density in the hadron, and we explain how to compute the fluctuations of the latter in the case in which the hadron is a heavy onium. In Sec. 3,we propose physical pictures of the events populating the tails of the gluon number distributions, and we establish analytical formulas to describe them. Section 4 contains the main new results of this paper, namely, a thorough numerical investigation of the tails of gluon density fluctuations in an onium. The final section 5 presents our conclusions. Technical details on the numerical calculations are gathered in the Appendices.

Multiplicity in hadron-nucleus collisions
We recall the picture of particle production in p-A scattering introduced in Refs. [21,20]. We first relate the final state particle multiplicity to the gluon number density in the Fock state of the incoming hadron at the time of its interaction with the nucleus, before explaining how the event-by-event fluctuations of the gluon density can be thought of.

Relation to the gluon number density in the hadron
Let us consider most generally the scattering of a dilute hadron, such as a proton or a quarkonium (which may be either a model for a hadron, or an actual state of a virtual photon), off a large nucleus, occurring at an energy corresponding to the total relative rapidity Y , assumed large compared to 1. The gluons in the Fock state of the hadron at the time of the interaction 1 that have a transverse momentum smaller than the saturation scale of the nucleus undergo scatterings, which may put them on-shell with high probability. The ones that have a transverse momentum larger hardly interact, and thus do not pick up the energy that would be needed to produce them: Therefore, they must recombine with other partons before they reach the final state. Hence, naively, the number of hadrons measured in the final state at a given rapidity y 0 with respect to the nucleus, in a given event, is proportional to the number of gluons with transverse momentum smaller than the saturation momentum Q s (y 0 ) of the nucleus in the corresponding Fock state of the hadron [20]. It is tantamount to the gluon number density integrated up to this momentum which carry a specific momentum fraction of the hadron in the initial state. We shall denote it by xG(x, Q 2 s (y 0 )). The ordinary gluon density xG would be equal to the mean of xG when averaged over the events.
More precisely, let us call M the mass of the onium and dN/dy the number of gluons per unit rapidity observed at an angle corresponding to the rapidity y 0 relative to the nucleus (resp. y ≡ Y − y 0 relative to the onium). Then, a calculation in the double-logarithmic approximation of QCD leads to [23,20] A schematic representation of the mechanism behind the correspondence formalized by this equation is given in Fig. 1. The gluons/hadrons produced in the final state have transverse momentum of the order of Q s (y 0 ): Indeed, the multiple scatterings broaden the transverse momentum of the gluons to that value. Note that in the present model, the saturation scale Q s (y 0 ) is a momentum which fully characterizes the nucleus [24] and depends only on the rapidity y 0 , not on the considered event. Indeed, since the nucleus in its ground state is already a dense object, the statistical fluctuations of its partonic content can be neglected throughout its rapidity evolution. The Fock state of the hadron instead results of a stochastic evolution up to the rapidity y, starting with a few partons. Hence it shows large event-by-event fluctuations, which directly translate into large fluctuations of the multiplicity observed in the final state.
Although our paper has a purely theoretical scope, a comment on the phenomenological applicability of our results and a justification of the relevance of our picture for proton-nucleus scattering at the LHC is in order. Our calculation requires a large nuclear saturation momentum Q s (y 0 ), of a few GeV, in order for perturbation theory to be justified; Hence the rapidity y 0 Initial state Scattering with the nucleus Final state Figure 1: Space-time picture of a scattering event of an onium off a nucleus. The frame is chosen such that the nucleus is leftmoving and boosted to the rapidity y 0 , and the onium is rightmoving at rapidity y ≡ Y − y 0 . In the particular event represented here, the onium fluctuates into 4 dipoles. Two of them are much larger than 1/Q s (y 0 ) and interact with the nucleus, while the two others are much smaller than 1/Q s (y 0 ) and do not interact. The scatterings are represented by gluon exchanges with the bulk of the nucleus. In a double-logarithmic scheme, only the gluons which are at the endpoints of the dipoles that scatter materialize in the final state, eventually in the form of hadrons of transverse momentum of the order of Q s (y 0 ). The other gluons recombine before they reach the final state. (The decay products of the nucleus are not represented).
should be relatively large. The hadron carries the remaining available rapidity, Y − y 0 , which is assumed to be large, but small compared to rapidities at which nonlinear effects should be taken into account in its evolution (namely, parametrically, Y − y 0 ln 2 1/α 2 s , see for example [25,26,27]). This is consistent with the kinematics of LHC, where the beam rapidity is of order Y 10, in such a way that we may pick a suitable value of y 0 , to allow a large-enough saturation momentum in the nucleus, and still moderate evolution in the proton.
The simplest hadron we can think of is a quark-antiquark pair in a color singlet state, which we call "onium". The Fock states of such an object are most easily described, analytically and numerically, in the framework of the color dipole model. Let us recall briefly how this works.

Event-by-event fluctuations of the gluon number density
When probed in its restframe with a wave of wavelength of the order of its spatial extension, an onium is viewed as a bare quark-antiquark pair: Indeed, although they are ubiquitous, the quantum fluctuations are too short-lived on the scale of the interaction time to play a role in the interaction. If instead the onium is probed with the same wave in a frame in which it has a large rapidity, then the lifetimes of its quantum fluctuations are Lorentz-dilated: Therefore, it appears essentially as a set of a large number of gluons.
One can evaluate the probability of a particular Fock state at a given rapidity by computing all diagrams contributing to the probability amplitude of finding the onium in that state. In the limit of a large number of colors N c , and in the leading logarithmic approximation (LLA) in which one keeps only the contributions for which the number of powers of y accompanying each power of α s is maximum, a convenient way to organize the calculation is the so-called color dipole model [28].
The dipole model uses coordinates in the two-dimensional plane orthogonal to the worldline of the onium, instead of momenta, to label the partons in the Fock state. Thanks to the large-N c limit, the set formed by the initial quark-antiquark pair along with its gluon fluctuations can be replaced by a set of color dipoles, the endpoints of which coincide with the position of the quark, of the antiquark, or of one of the gluons. The graphs contributing to a given state are generated by a stochastic branching process in rapidity [28]. The latter is completely defined by the elementary probability that a dipole defined by the pair of the two-dimensional position vectors of its endpoints, (x 0 , x 1 ), branches into two dipoles, (x 0 , x 2 ) and (x 2 , x 1 ) respectively, by emitting a gluon at position x 2 when its rapidity is increased by the infinitesimal amount dy. A calculation in the framework of perturbative QCD leads to the following expression for the probability [28] dp 0 | pQCD =ᾱdy whereᾱ ≡ α s N c /π, and we introduce the notation Analyzing this equation, one sees that there is a non-negligible probability that the gluon be emitted at a large distance of the initial onium. However, confinement should forbid, at least in principle, the production of dipoles which are bigger than, typically, 1/Λ QCD . Being an intrinsically nonperturbative effect, we cannot attain it through a perturbative calculation. Therefore, we add it to the original model in the form of a "cutoff function" Θ that forbids dipoles of size typically larger than some infrared length scale R ∼ 1/Λ QCD to be produced. This leads to the following modification of the splitting probability: Most generally, the cutoff function Θ must satisfy the following limits: It is also expected to reach 0 "fast enough" (that is, at least exponentially) when x 02 or x 12 become larger than R. This function is arbitrary in our treatment, however, it will become clear that the observables we are interested in can depend only marginally on its precise form. The relation between the number of dipoles and the gluon density is very simple if one restricts oneself to double-logarithmic accuracy: where n(r s ; x 01 , y) is the number of dipoles of size larger than r s in the state of an onium of initial size x 01 , observed at rapidity y. The derivative enters the right-hand side because xG(x, Q 2 s ) is the density of gluons of a fixed momentum fraction x, while n(r s ; x 01 , y) enumerates the dipoles which have a rapidity smaller than y. It is the double-logarithmic approximation that enables one to identify the size r s to the inverse momentum 1/Q s ; Sizes and momenta being conjugate to each other through Fourier transform, this identification does of course not hold in general.
Thanks to Eq. (5), xG, and thus, through Eq. (1), the number of particles produced in a given rapidity slice in the final state, have the same fluctuations as the number of dipoles in the Fock state of the onium at the time of the interaction. Therefore, the scope of the following sections will be to study first analytically, and then numerically, the probability P n (r s ; x 01 , y) to have n dipoles of size larger than r s in the Fock state of the onium after evolution of a dipole of initial size x 01 over the rapidity interval y.

Tails of the dipole number distribution
In this section, we study analytically the high and low-multiplicity tails of the dipole number distribution P n (r s ; x 01 , y), developing physical pictures which will prove useful for the interpretation of the numerical data.
"High" and "low" are intended with respect to the expected multiplicity. In both cases, we will assume that the dipole numbers are much larger than unity. In this limit n 1, the probability P n , which is defined as a function of the integer n, can be thought of as a function of a continuous variable, and thus as a probability density. The probability to observe a number n of dipoles in the interval [n 1 , n 2 ] then reads n 2 n 1 dn P n .

Heuristics
No infrared cutoff. Let us recall that in perturbation theory, in the absence of an infrared cutoff, the rapidity-evolution of the expected number of dipoles larger than some size r s , starting from an onium of size x 01 , is governed by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [29,30] ∂ ∂y n (1) (r s ; x 01 , y) = dp 0 | pQCD dy n (1) (r s ; x 02 , y) + n (1) (r s ; x 12 , y) − n (1) (r s ; x 01 , y) , (6) where the integration goes over the whole transverse plane. Denoting byᾱχ(γ), where the eigenvalue of the kernel of the BFKL equation associated to the eigenfunction (x 2 01 /r 2 s ) γ , we can write the solution of equation (6) as a continuous superposition of the eigenfunctions weighted by eᾱ yχ(γ) . The initial condition corresponding to one single dipole of size x 01 reads n (1) (r s ; x 01 , y = 0) = θ(x 01 − r s ). The solution to the BFKL equation then reads The leading behavior of n (1) at large rapidities is given by a saddle point: It is useful to represent the dipole evolution in the two-dimensional (r,ỹ) plane as the curve connecting the points (x 01 , 0) and (r s , y) and that solves the saddle-point equatioñ when the value of γ s is fixed by the boundary conditions, i.e. when it solves Eq. (9). In log r scale, this curve is just a straight line, see Fig. 2. It represents the most probable evolution path.
If instead of the expected dipole number n (1) we are interested in the probability P n of having a given number n of dipoles in the Fock state, or if, equivalently, we focus on the set of moments n (k) of the dipole number, then the path that corresponds to the main contribution is not necessarily a straight line. It was shown in Ref. [21] that when n is large compared to its expected value (or equivalently, k is large compared to 1), this path essentially consists in two steps: The initial onium generates, through a fast evolution, a dipole of large size r max , which subsequently decays into many (mainly smaller) dipoles. The presence of the large dipole at an early stage of the evolution is necessary if one asks for a large multiplicity, because large dipoles yield much more offspring of size larger than r s than smaller ones. This first step has a low probability, which decreases as r max increases. But on the other hand, the number of offspring increases with r max . The optimal size r max of the intermediate large dipole depends on the maximum rapidity and of the final number of dipoles of which the probability is evaluated. Following Salam in Ref. [31], we assume that the production of the large dipole occurs literally in the very first step of the evolution. As we will check a posteriori, its size increases with n, and thus r max can be made arbitrarily large, say r max r s , by selecting very large values of n at fixed y. Once the large dipole has been produced, it decays into much smaller ones. This second step in the evolution is dominated by decays which are strongly ordered in the dipole sizes, from large to small. The solution to the saddle-point equation in Eq. (9) is close to γ s 0, a region in which χ(γ s ) may be approximated by 1/γ s . In this limit, the number of dipoles larger than r s resulting from the decay of a dipole of size r max r s reads This is the well-known "double-logarithmic" limit. In Ref. [21], it was proven that the fluctuations of the dipole number n on such an evolution path are suppressed exponentially. We will check a posteriori that these fluctuations are overall negligible, and thus, that we can assume that the second step of the evolution is deterministic. Consequently, the probability to observe more than say N dipoles, ∞ N dn P n , coincides with the probability that r max be larger than R, where R is such that n (1) | DL (r s ; R, y) = N . Solving this elementary equation for R, the relation between the probability distribution of the dipole number n to that of the size of the intermediate large dipole takes the following form: The probability that the initial dipole of size x 01 splits into a dipole of size r larger than some given R, itself much larger than x 01 , is suppressed by the ratio of the squared sizes, see Eq. (3). Indeed, the distribution of the sizes of dipoles produced in the splitting of a dipole of size x 01 , conditioned to the occurrence of a splitting into similar or larger-size dipoles reads 1 Nᾱ dp 0 | pQCD dy , where N is a normalization factor of order 1. Thus the probability of having a dipole of size r max larger than R is just the following integral: Now, to arrive at the distribution of the dipole number n, it is enough to replace R by r s e ln 2 N/(8ᾱy) (see Eq. (12)) in the previous equation. Note that this quantity is the typical size of the intermediate dipole, which confirms the a priori assumption made above that it grows with n. Taking finally the derivative with respect to N and evaluating it at N = n, we obtain the density 3 P n : We note that the large-n tail of P n is much fatter than an exponential decay. This is the a posteriori justification for having neglected the stochasticity in the second step of the evolution, consisting in the decay of the large dipole of size r max . A comment is in order. Dipole evolution is often assimilated with a branching random walk (BRW). This is correct when one looks, for example, at an observable probing the number of dipoles overlapping with a given point in the transverse plane. But in a BRW, the number of objects after some given evolution is distributed exponentially, by contrast with Eq. (14). The fat tail we find here is a feature of QCD which shows up when we count all dipoles (larger than a given size, to talk of an infrared-safe quantity) independently of their transverse position. Technically, it is related to the size-dependence of the dipole splitting rate, while in a BRW such as e.g. the branching Brownian motion, particle splitting and diffusion are completely uncorrelated. We are going to see that we actually recover an exponential distribution when we enforce an infrared cutoff on the evolution.
Enforcing an infrared cutoff. We now consider the modified dipole model which incorporates an infrared cutoff in the form of the Θ function, see Eq. (3).
With an IR cutoff, the typical size of the dipoles generated in the first step of the evolution eventually becomes limited by the infrared boundary if one focuses on very large values of n. Hence the size r max of the intermediate dipole will always end up being of order R. So unlike in the purely perturbative QCD case, the stochasticity in n cannot come from the first step. It necessarily stems from the second step, consisting in the decay of the large dipole. Let us try to understand the form of the distribution of these fluctuations.
The IR cutoff forces the produced dipoles to be smaller than, typically, R throughout the evolution. On the other hand, because it is probabilistically disfavored, a small dipole does not split to much larger dipoles. So starting from a dipole of size close to R, the final number of dipoles is essentially built up by a backbone of successive splittings of dipoles to similar-size or smaller, but not much smaller, dipoles, each of which gets dressed by a number of very small dipoles (of size of order r s ) proportional to the rapidity interval over which it evolves. Hence this second step in the evolution essentially looks like a 1 → 2 branching process, in which the branchings occur at an almost constant rate, as in a BRW (see the comment above). The fluctuations in such a process are known to be exponential, ∝ e −n/n 1 /n 1 , where n 1 is the mean number of objects eventually produced at the final rapidity. 4 Hence, we expect the shape of the dipole number distribution to follow such a law. We now need to understand the overall normalization, as well as the parameter n 1 .
Concerning the slope of the exponential, n 1 , a good estimate can be obtained from the mean number of dipoles produced by an initial dipole of size r max of order R. It satisfies a modified BFKL equation, i.e. Eq. (6) with the substitution dp 0 | pQCD → dp 0 , which however cannot be solved exactly because the eigenfunctions of the kernel of such an equation are not simple in general. However, we may obtain a good approximation to its solution by replacing the Θ function by a sharp Heaviside distribution, which in turn is tantamount to an absorptive boundary. Then, the method of images can be used to arrive at a solution to this problem.
We start with the solution to the ordinary BFKL equation without a cutoff, Eq. (8). We may evaluate the integral in the saddle point approximation for large y, using Eq. (9) with the substitution x 01 → r. We anticipate that the saddle point equation (9) for γ s has a solution near γ s = 1 2 , therefore we replace χ(γ) by its expansion around γ = 1 2 : where ζ is the Riemann zeta function. Then Following Ref. [33], the absorptive boundary is implemented through the method of images applied to the diffusive part, represented (up to numerical constants) by the factor in the curly brackets in the previous equation. The result reads . (17) It is easy to check that this solution obeys the (ordinary) BFKL equation, and that it satisfies indeed the boundary condition n (1) | sp,Θ (r s ; R, y) = 0. We observe that there is an optimal dipole size, that maximizes the mean number of dipoles at the end of the evolution: n (1) | sp,Θ as a function of r max exhibits a maximum at r max = O(R), located between 0 and R. Indeed, it vanishes linearly with r max as r max → R as a consequence of the presence of the absorptive boundary, and also goes to zero as r max → 0.
Concerning the normalization of the exponential, the probability to generate a dipole of size r max of order R in the first step of the evolution is a factor in P n independent of n for large n. It can be estimated by replacing R by R in Eq. (13): The result is proportional to x 2 01 /R 2 . All in all, these heuristic considerations lead us to the following expression for P n : This formula is indeed of the same form as the one derived in Ref. [21] with a different, more mathematical, method. The advantage of the present heuristic approach is that it comes with a simple picture of the evolution of the Fock states into high-multiplicities, while the more abstract method uses the factorial moments, for which it may be more difficult to build an intuition. Also, we see that we have not used any detailed property of the cutoff function Θ: It just needs to be "sharp enough", namely it must decay faster than some power of R/x 02 , when x 02 (and thus x 12 ) gets larger than R. This shows that the high-multiplicity tail of the distribution of n cannot be very sensitive to the precise form of Θ.
Note that there is actually an awkward point in our heuristic discussion. Indeed, when the initial dipole has a small size x 01 compared to the infrared cutoff R, if the production of the large-size dipole really consisted in one single splitting, then due to the geometry of dipole splitting, one would have two dipoles of sizes x 02 and x 12 very close to each other, and not one. 5 Indeed, x 01 = x 02 + x 21 , so if x 02 ∼ R x 01 , then x 12 ∼ x 02 . But in a situation in which x 02 = x 12 , a simple calculation shows that the fluctuations of the number of offspring of this pair of dipoles would be distributed as ∝ n × e −n/n 1 instead of a simple exponential. As we will see in the detailed simulation of Sec. 4, this would contradict our numerical results. There may be two ways out. First, the two dipoles are never exactly of the same size, and consequently, the mean dipole yields n 1 associated to each of them are not exactly the same (see Eq. (17)). Then, for very large n, the fluctuations are always dominated by the offspring of one of the dipoles (the one that yields most offspring on the average). Second, in the more detailed analysis of Ref. [21], the production of the large dipole resulted from a BFKL-like evolution, not from one single splitting (although that evolution turned out to be very fast when n was set to be very large). In this case, there is no reason why there should systematically be two large dipoles of (almost) identical size. Finally, the exponential decay (without a n-dependent prefactor) of P n with n was found in the more straightforward calculation presented in Ref. [21] (and reproduced, for completeness, in the next section), and it seems well supported by the numerical data, see below Sec. 4.

Solution from an Ansatz.
We introduce the generating function Z(r s ; x 01 , y|u) of the factorial moments of the dipole number: It is well-known that it obeys the Balitsky-Kovchegov (BK) equation [14,36] (modified by the infrared cutoff here): Θ(x 02 , x 12 ; R) [Z(r s ; x 02 , y|u)Z(r s ; x 12 , y|u) − Z(r s ; x 01 , y|u)] .
(20) For the analytic calculation, we choose a factorized form for Θ: where the functionθ has the limits θ(X) Its precise form is not really relevant for the asymptotic calculations we will carry out, except for one step (see below), for which we will need to pick a specific function forθ, to arrive at a simple expression. Even in the case of purely perturbative QCD (R = +∞ or equivalentlyθ = 1), we do not know how to solve the BK equation (20) accurately enough to be able to extract the probabilities P n from the solution for Z. However, we may notice that the large-n asymptotics of P n is connected to the large-k asymptotics of the factorial moments n (k) . It is straightfoward to convert the BK equation into a hierarchy of equations for n (k) : where the factorial moments n (k) are defined as the event-averages of the products n(n−1) · · · (n− k + 1), i.e.
Again, it is not possible to find an explicit solution for n (k) -unsurprisingly, since the infinite set of equations (23) is equivalent to Eq. (20). However, it is not difficult to figure out a plausible Ansatz. We try where the C k 's are constants. Taking n (1) = n (1) | sp from the saddle-point solution of the BFKL equation, Eq. (9), inserting Eq. (25) into Eq. (23), and keeping only the leading term in the limit of large rapidities, n (1) ∝ eᾱ yχ(γs) , the equation for the moments n (k) boils down to an equation for the constants C k : The integrals over x 2 are finite functions of x 01 /R. We anticipate that the first terms in the right-hand side, inside the curly brackets, are negligible compared to the other terms. Under this assumption, which we will check a posteriori, the equation to solve simplifies to To push further the analytical calculation, we now need an explicit form for the infrared cutoff functionθ. We choose a Gaussian:θ (X) = e −X 2 /2 .
After the appropriate replacements have been done, the integration over x 2 in Eq. (27) can be performed: which is justθ(x 01 /R)/4. Therefore, Eq. (27) boils down to an algebraic recursion for the constants C k of the form Inspection of this equation straightforwardly shows that for asymptotically large k, C k behaves like a k k!, where a is a constant.
We are now in a position to go back to Eq. (26) and check a posteriori that it was indeed justified to neglect the terms in the curly brackets. This simply stems from the fact that these terms are overall proportional to C k , while the sum of all the other terms, the ones we have kept, gives a contribution proportional to k × C k , which is much larger for k 1. Putting everything together, we see that the factorial moments of the dipole number read where n 1 = a × n (1) and c is a constant which we expect to be of order 1. Now, from the knowledge of the moments we can obtain the probability density function P n . Since we deal with typical multiplicities much larger than 1, factorial moments of order k can be approximated by ordinary moments of the same order, i.e., which, using Eq. (31), leads to the final result: This is fully consistent with Eq. to split into large dipoles until a relatively large rapidity is reached. The forbidden region is represented by the shaded area. In that region, the evolution is limited to the emissions of small dipoles, of size smaller than r 1 .

Low-multiplicity tail [37]
We now turn to the evaluation of the distribution of the low multiplicities. By "low" we mean much lower than the typical or mean multiplicityn n (1) , but at the same time still much larger than 1. This region has not drawn as much attention as the high-multiplicity region, with the exception of Ref. [38], in which Iancu and Mueller analyzed it with a view to understanding the Levin-Tuchin law [39] for total dipole-nucleus versus dipole-dipole cross section deep in the saturation region.
The only way to generate events with n n is to veto the splittings of the initial dipole in the beginning of the evolution, except if the latter are small enough: Indeed, we know that dipoles which have sizes close to the saturation radius r s cannot evolve into high-multiplicity states except by creating large dipoles, but this has a large cost in probability which makes such a process subdominant. On the other hand, once the initial dipole has split into similar-size or larger dipoles, then the cost of keeping the density of the state low becomes large. Once a couple of dipoles have been emitted, the subsequent evolution can be considered deterministic, and the latter generates a number of dipoles which grows fast with the rapidity. Hence we expect the low-multiplicity tail of P n to be made of events in which the occupation number is kept low throughout the initial stages of the evolution. The evolution of such configurations is schematized in Fig. 3.
We are going to derive an expression for the low-multiplicity asymptotics of P n from these simple considerations. We shall assume that the probability of a given dipole number n coincides with a suppression factor for dipole splitting inside an appropriate region D of rapidity and transverse size, up to slowly-varying prefactors that we shall discard: For n small enough compared to the typical multiplicityn, D must also include relatively small dipoles compared to x 01 , namely either x 02 x 01 or x 12 x 01 . Let us call r(ỹ) the lower boundary of D at some fixedỹ, namely the minium size of the dipoles included in the domain over which we integrate. Since the integral over x 2 diverges logarithmically when r(ỹ) goes from O(x 01 ) to zero, while the contribution of the dipoles larger than x 01 to the integral is finite and of order 1, keeping only the strongly-ordered regions x 02 x 01 , x 12 ∼ x 01 and x 12 x 01 , x 02 ∼ x 01 is enough to get the dominant term in the limit r(ỹ) x 01 . We write Coming back to Eq. (34), changing variable from x 02 to ρ = ln x 2 01 /x 2 02 , the following approximation can be written for the probability: We pick the simplest for D: We assume it just consists in the rectangular-shaped region (r,ỹ) ∈ [r 1 , +∞[×[0, y 1 ]. We choose y 1 and r 1 such that a dipole of initial size x 01 starting to evolve deterministically at y 1 produces exactly n dipoles at the final rapidity y, and a dipole of size r 1 starting to evolve at rapidity 0 also produces n dipoles at y. It is not difficult to figure out that these conditions are enough to guarantee that no dipole emitted outside of D can grow into a state of multiplicity much larger than n. Hence in the presence of such a vetoed region, writing ρ 1 = ln x 2 01 /r 2 1 , the distribution of the number of particles reads P n ∝ e −ᾱy 1 ρ 1 .
Note that by choosing a rectangular region D, we neglect a term in ln P n which is proportional to (ᾱy 1 ) 2 . We now use the defining conditions for y 1 and r 1 to express these variables with the help of n andᾱy. In the double-logarithmic approximation (11), these conditions read e 2 √ᾱ (y−y 1 )ρs = n and e 2 √ᾱ y(ρs−ρ 1 ) = n, where we introduced the notation ρ s = ln x 2 01 /r 2 s . The previous equations enable us to rewrite Eq. (37) as Finally, for the sake of writing down a more compact formula, we may use again the doublelogarithmic approximation to express the productᾱyρ s with the help of the expected dipole numbern: 2 √ᾱ yρ s lnn. The following expression is obtained: or, when expanded and ordered by decreasing importance in the limit n n: ln P n = 1 2 ln 2 n − 1 4 A remarkable feature of this result is that the leading n-dependence at fixed n and large rapidities does not involve at all the infrared, and actually does not depend on any scale at all.

Numerical study
In this section, we test the validity of the physical pictures proposed in Sec. 3 to understand the behavior of the tails of the multiplicity distribution. To this aim, we perform high-statistics Monte Carlo simulations of dipole evolution for different values of the parameters. We measure distributions of the dipole multiplicity, P n , and compare their large and the small-multiplicity tails to Eqs. (33) and (41) respectively, for different values of the parameters. Although we do not report on it here, we have also tested many choices of a rapidly falling Θ, and checked that the qualitative shape of the tails were not altered [40], as expected from general considerations. In the numerical results we will present, we will restrict ourselves to the Gaussian IR cutoff (Eq. (28)) which was employed in the analytical calculations of Sec. 3.

High-multiplicity tail
The parameters of our numerical simulations are set to be the following: r s /R = 1/40, x 01 /R are varied between 0.67 and 5 × 10 −2 , andᾱy = 2 to 5. Our analytical results rely on the fact that the mean evolution between R and r s is driven by an eigenvalue of the BFKL equation close to χ( 1 2 ): In other words, the solution to the saddle point equation in Eq. (9) (with x 01 → R) was assumed to sit around γ s = 1 2 . Let us check that it is indeed the case with the set of parameters we have chosen: (42) whenᾱy is set to be the value for which we have collected most of the data, namelyᾱy = 4.

Higher-order moments
To study the behavior of the high-multiplicity tail, we look at the higher-order moments of P n . In particular, for a given moment n (k) , we look at its behavior as function of k, for different values of the size of the initial dipole, x 01 . Our goal is to check that P n has an exponential tail, and to provide a measurement of its slope. If P n were a strict exponential distribution of the form, say, P exp n = e −n/n 1 /n 1 , then its moments would simply read the ratio of successive moments would satisfy the following equality: Since the large-multiplicity tail of the distribution is probed by moments of high order, we look at the behavior of this ratio for large k.
It is instructive to estimate the typical value of n s probed by a moment of order k. This is given by the value of the dipole number, n, that contributes most to the integral in Eq. (43). When k is large, a saddle point at n s = kn 1 dominates the integrand. These are the typical values of n probed by the k-th moment.
Numerical results for Eq. (44) in the case of initial dipoles of different relative sizes x 01 /R, evolved up toᾱy = 4, are shown in Fig. 4. The very smooth behavior of the data points and of the magnitude of the statistical errors for different values of k is due to an intrinsic correlation of errors at different k, which comes from the fact that for a given value of x 01 /R, n (k) 's are computed using the same sample of data. We draw a horizontal line at n 1 = 10500 as an illustrative value of n 1 that is compatible, within one sigma, with all the curves shown in the plot for large-enough k. We conclude that this value is independent of x 01 /R.
We now check that the rapidity dependence of n 1 is essentially exponential, up to prefactors, as predicted by Eq. (9). For selected values of x 01 /R, we report results at different values ofᾱy in Fig. 5. Note that for any value of y, this ratio tends to a constant n 1 (y), which is independent of x 01 /R. The logarithmic scale on the y-axis makes it obvious that n 1 grows with y approximatively like an exponential, in agreement with the asymptotic identification n 1 ∼ n (1) .
Finally, our Ansatz predicts that the coefficient multiplying the exponential in Eq. (33) should present a specific quadratic dependence on the size of the initial dipole, x 01 , for x 01 R. In order to test this, we exploit the fact that this coefficient appears in the expression of the moments, Eq. (25). This implies that the ratio of two moments n (k) computed at two different values of x 01 /R provides direct information about this coefficient. Hence, we compute the ratio n (k) [x 01 /R]/n (k) [x 01 /R = 1/6] in our calculations atᾱy = 4. The results are shown in Fig. 6 up to k = 12, after which statistical uncertainties dominate. We compare the numerical data with both a simple quadratic behavior (dashed line), and a quadratic Ansatz corrected with a Gaussian factor (solid line), i.e., the full prefactor in Eq. (33). We observe that the data points tend to fall on the expected curves 6 as we move to larger values of k. Moreover, it is clear that the dashed and the solid line describe equally well the trend of the data points in the region where x 01 /R is not close to unity. We emphasize that this result is very important, because it implies that the simple heuristic discussion leading to Eq. (18) allows us to correctly predict   the behavior of the numerical calculations. This confirms the robustness of our intuitive picture of high-multiplicity events, and supports our statement that the exponential tail is universal irrespective of variations of the (sharp) infrared cutoff function.
In summary, all the expectations of Sec. 3 are confirmed by our numerical results. We conclude that, within our uncertainties, the high-multiplicity tail is an exponential correctly captured by Eq. (33).

Shape of P n at large n
Let us show, then, the actual shape of the distributions of dipole (gluon) number obtained in our Monte Carlo calculations. Results are shown as circles in Figure 7. We exploit the result obtained in the previous subsection to characterize the tails at high multiplicity. We overlay our curves with the asymptotics shown in Eq. (33), where we set γ s = 1 2 , and use n 1 = 10500, i.e., the value inferred from the analysis of Fig. 4, and we choose the overall normalization factor to be c = 2. The asymptotic curves are shown as red dashed lines in Fig. 7. We find that, with common values for n 1 and for the normalization, c, the exponential asymptotic trend is able to describe all the multiplicity distributions at large n, irrespective of x 01 /R. We note that the constant c is indeed of order unity, as expected from the theory.  numerical data for n 1 . Putting numbers into Eq. (17), we find that forᾱy = 4 and R/r s = 40, n 1 reaches a maximum of about 11700 for r max 0.37 × R. This is consistent with the value of n 1 found in the numerical data (n 1 ∼ 10500). Remarkably, it is for such values of x 01 that P n is closest to a pure exponential, see Fig. 7. This means that when one sets x 01 to the value r max expected to maximize the number of final dipoles, then the optimal path in the (r,ỹ) plane (see Sec. 3) is a straight line. The emission of large dipoles in the initial steps is no longer advantageous, since the probability is strongly dumped by the cutoff function. This is a nice consistency check of the whole picture.
We have checked that the y-evolution of n 1 which we may deduce from Fig. 5 is also qualitatively reproduced by Eq. (17).

Low-multiplicity tail
Now we turn our attention to the low-multiplicity tail of the dipole number distribution, P n . Again, we want to check that the physical picture described in Sec. 3 for the low-n tail is consistent with the Monte Carlo results. To this aim, we shall perform fits of P n in the lowmultiplicity region using the formula in Eq. (41). Figure 8(a) shows the results for the low-multiplicity tails of the distributions obtained at different values of x 01 /R, forᾱy = 4. The dashed lines in the figure represent fits obtained using Eq. (41). The strategy of our fitting analysis is the following: For each tail, we fit the same number of points (around 15), starting from the points with lowest probability. We checked that, doing so, we do not break the condition n n, which is required for Eq. (41) to apply. We find that the quality of the fit improves as we move to larger values of x 01 /R, as the χ 2 per degree of freedom (dof) of the fit is close to 4 at x 01 /R = 0.05, and close to unity at x 01 /R = 0.5. This behavior confirms our expectation that the physical picture presented in Sec. 3 works as long as 1 n n, a requirement which is loosely satisfied by the curves at small x 01 /R, where we observe events with n ∼ 1. For each tail, the two-parameter fits return one overall normalization, and a value for the mean of the distribution,n. Remarkably, we find that the overall normalization varies by less than a factor 2 among the different fits, as long as x 01 /R > 0.1. This supports our conclusion of a universal shape for the low-multiplicity tail, which is simply shifted towards larger values of n ifn increases. As for the fitted values ofn, we obtain numbers which are smaller than the actual mean value of the distribution, although in a systematic way: The fittedn at x 01 /R = 0.5 is much closer to the true value than for the fit at x 01 /R = 0.05. Again, we understand this as a consequence of the fact that in our calculations we do never reach the fully asymptotic regime 1 n n, since even for x 01 /R = 0.5 we observe events with n ∼ 10. These systematics support the validity of Eq. (41).
To corroborate this statement in more generality, it is useful to look at the results shown in Fig. 8(b). In this panel, we focus on the multiplicity distribution at x 01 /R = 0.05, which is the distribution providing the least satisfactory results from the fit of the low-n tail. We show how the tail evolves with the rapidity,ᾱy. We see that increasing rapidity has an effect similar to that of increasing x 01 /R, in the sense that the distribution gets shifted towards larger values of n. Clearly, the quality of the fit improves as we increase the value ofᾱy, except for the largest αy = 5, for which the trend of the data is less well reproduced by the theoretical curve. This is actually expected, because the double-logarithmic approximation, crucial in the derivation of Eq. (41), is less justified for larger values ofᾱy.
We now show that the shape of P n in the low-multiplicity region is not affected by the presence of the infrared cutoff, Θ. To this aim, we perform calculations with fixed ratios x 01 /r s and x 01 /r UV , whereas we vary R/x 01 , i.e., the distance between the initial dipole size and the infrared cutoff. We test very different scenarios: We consider a range of values for R/x 01 , as well as the case where the infrared cutoff is absent, i.e., R/x 01 = +∞. We obtain the results shown in Fig. 9. The two-parameter fits, obtained from Eq. (41), are of excellent quality, and return a χ 2 /dof close to unity. We conclude that the infrared cutoff has simply the effect of translating P n towards lower values of n, without altering the shape of the low-multiplicity tail. Note that the agreement we have found is better than one would have expected, considering the values of the parameters we have chosen. Indeed, the double-logarithmic approximation was heavily used in deriving Eq. (41), and the latter is supposed to be best justified for γ s → 0, where γ s is evaluated using Eq. (42), with the substitution R → x 01 . It turns out that γ s ranges between 0.41 and 0.46 when computed with the parameters with which our data has been generated.
Therefore, it is important to confirm that the analytical formula (41) works also very well in the kinematical region which satisfies the assumptions made for its derivation. To this aim, we generate data at lower rapidities and much larger values of x 01 /r s : αy = 1 and x 01 r s = 10 8 =⇒ γ s 0.16 αy = 2 and x 01 r s = 10 5 =⇒ γ s 0.28.
The corresponding values of γ s are now such that the DL approximation is very well justified. The agreement between the numerical data and the theoretical parametrization is indeed extremely good. Both fits return a value of χ 2 /dof close to unity; see Fig. 10.
In conclusion, the shape of gluon number distribution in the region of low multiplicity is not affected by the presence of an infrared cutoff, and it is correctly captured by the physical picture of Sec. 3.2, as long as the parametersᾱy, x 01 and r s are such that there is a large-enough region in which the inequalities 1 n n are satisfied.

Summary and outlook
We have studied analytically and numerically the tails of the multiplicity distribution of gluons corresponding to dipoles larger than a given size r s ∼ 1/Q s (y 0 ) in the Fock state of an initial onium.
The low-multiplicity tail shows very robust features, as its leading dependence in n does not depend on any scale or parameter. In particular, it does not depend on the infrared cutoff, and we have been able to describe it with Eq. (41) even at the border of the double-logarithmic region. The large-multiplicity tail is exponential: A proxy for the parameter of this exponential is the mean integrated gluon density in a dilute hadron evaluated at the saturation scale Q s (y 0 ). The onset of the exponential depends on the ratio of the onium size to the confinement size.
It is important to appreciate that the fluctuations of the gluon density, which are relevant for the multiplicity distribution in hadron-nucleus scattering we have studied, are fundamentally different in nature from the event-by-event fluctuations of the saturation scale discussed in Ref [41,42]. The source of these fluctuations was identified to be the saturation of the dipole density in the evolution itself (whose mechanism may be recombinations, or screening of further emissions), which we have not implemented in the dipole model and which are not relevant at the energies of the present colliders. The fluctuations studied here are also different from the so-called "front" and "tip" fluctuations in branching random walks identified in Ref. [43,35]. The latter prove relevant in QCD for observables for which only dipoles overlaping with a given impact parameter may play a role, or only the fluctuations of the size of the largest dipole(s): The total and diffractive γ * -A deep-inelastic scattering cross sections are examples of such observables, whose evolution can be considered in the universality class of branching random walks. The fluctuations which are at work in the process we have been analyzing in this paper are essentially equivalent to the statistical noise of the total number of particles generated by a 1 → 2 branching process after some given evolution.
Going from our study of the multiplicity of gluons produced in onium-nucleus collisions to the experimental data on the hadron multiplicities in proton-nucleus scattering will require some modeling. First, the details of the state of the initial proton will certainly determine the overall normalization, but may also introduce n-dependent prefactors. Second, a hadronization model will be needed to link our parton-level calculation to the actual hadronic final state. The modelization of the proton is presumably not easy. Several recent studies model its color field as a classical Gaussian field, see Ref. [44] and references therein. We deem that a set of a few partons is more appropriate, but whether a diquark dipole would be a good representation for its groundstate is questionable. Exclusive diffractive processes may help constraining the initial condition for the evolution, see e.g. Ref. [45].
But while the detailed form of the multiplicity distribution will certainly depend on the model for the initial hadron, we expect its general features to be robust. An interesting observable that could be amenable to dipole model description is multiplicity in the final state of deep-inelastic γ * -A scattering, for sets of events in which the photon virtuality is chosen smaller than the nuclear saturation scale. In such a case, the virtual photon interacts through its quark-antiquark pair fluctuations, the distribution of the size of which is obtained from a straightforward QED calculation. Therefore, this process is very close to onium-nucleus scattering studied in the present paper.

A Implementation of the modified dipole model
In this appendix, for completeness, we describe our numerical implementation of the color dipole model. It actually follows closely the one in Ref. [31]. Other versions of the code have been written, incorporating different variations of the model to make it more suitable for phenomenological studies, such as energy conservation or gluon recombination; see e.g. Ref. [46].
The dipole model is a 1 → 2 branching process. Each dipole of a given set may split, independently of the other dipoles, into two dipoles when the rapidity is increased by dy. For definiteness, let us consider a generic dipole whose endpoints are labeled by the two-dimensional vectors x 0 and x 1 . The probability dp 0 it emits a gluon at position x 2 up to d 2 x 2 was given in Eq. (3): This probability diverges when x 2 coincides with the endpoints x 0 or x 1 : as well-known, the probability to split into very small dipoles can be arbitrarily large due to the collinear singularity. We thus need to introduce a lower cutoff r UV on the sizes of the produced dipoles in order to get a meaningful distribution. We choose to enforce it as sharp Heaviside θ function: where r UV is an arbitrary ultraviolet regulator, that needs to be taken much smaller than all distance scales relevant to our problem in such a way that it does not affect the physical results. The value of r UV we chose in practice is checked to satisfy this requirement in Appendix B.
Having the expression of the probability that the dipole splits in the infinitesimal rapidity interval dy, we can easily write the expression of the probability that the dipole splits (for the first time) at finite rapidity y (up to dy) by emitting a gluon at position x 2 (up to d 2 x 2 ): dp = dp 0 θ(x 02 − r UV )θ(x 12 − r UV ) e −ᾱλ(x 01 ;r UV ,R)y whereᾱλ is the inverse "lifetime" of the dipole, namely the inverse of the typical rapidity interval between two successive dipole splittings: λ(x 01 ; r UV , R) = 1 α dp 0 dy = d 2 x 2 2π x 2 01 x 2 02 x 2 12 θ(x 02 − r UV )θ(x 12 − r UV )Θ(x 02 , x 12 ; R).
We use two elementary techniques to implement dipole splitting in a Monte Carlo code. The first one is based on the following mathematics: If f (x) is the probability density of the real variable x, if F (x) = x −∞ dx f (x ) is its cumulative distribution function, and if y is distributed uniformly between 0 and 1, then x = F −1 (y) is distributed according to f . The algorithm that follows from this observation is practical whenever F and its inverse have analytical expressions. When this is not the case, then we use a rejection algorithm: We pick a densityf (x) whose inverse cumulative distribution may be expressed by a simple analytical formula, and which is such thatf (x) ≥ f (x) for all x. We then generate realizations of x according tof , and accept the generated values with probability f (x)/f (x).

A.1 Dipole evolution without an infrared cutoff
We first address dipole evolution defined by the probability dp in which the infrared cutoff is put to R = +∞, namely with the cutoff function set to Θ = 1: we shall denote this probability by dp| pQCD . We start by explaining how to generate the distribution of the position of the gluon (or, equivalently, of the size vectors of the produced dipoles in the splitting). Then, we generate the rapidity at which the splitting occurs.
Thus this law is completely determined by λ. Once this parameter has been computed, realizations of ∆y are trivial to generate. λ is given by a two-dimensional integral, over x 02 and ϕ: λ(1;r UV , ∞) = which needs to be performed numerically.
Generation of a full Fock state of an onium at rapidity y. Once one gluon emission, corresponding to a 1 → 2 dipole splitting, is implemented, the full Fock state at a given rapidity y is generated through a simple iteration. Starting from one dipole (x 0 , x 1 ), the rapidity ∆y at which it emits a gluon is generated. If this rapidity is found to be larger than the final rapidity y, then the evolution stops: the Fock state in this particular event consists in a single dipole. If instead it is less than y, then the position x 2 of the gluon is generated (see above), and the initial dipole is replaced by two dipoles, (x 0 , x 2 ) and (x 2 , x 1 ). This procedure is then just applied recursively to the two new dipoles over the rapidity interval y − ∆y.

A.2 Enforcing the infrared cutoff
When the cutoff R is set to a finite value, i.e. when a function Θ ≤ 1 is added to the splitting probability, then the integral over the angle ϕ can in general not be done analytically.
Two nested numerical integrals have now to be performed to compute the inverse lifetime λ. In practice, it proves useful to construct a lookup table containing the numerical evaluations of the two-dimensional integral for a set of sizes x 01 .
As for the distribution of the position of the emitted gluon, the simplest is to generate the dipole sizes with the weight given by the dipole model without the cutoff, and then use a rejection algorithm to correct the distribution: The splitting of a dipole of size x 01 into two dipoles of respective sizes x 02 and x 12 is accepted with probability Θ(x 02 , x 12 ; R). Since the cutoff function Θ cuts out a low-probability region, in which two dipoles larger than the initial one are produced, the efficiency of such an algorithm is quite high.

B Numerical test of the effect of the UV cutoff
The ultraviolet cutoff protecting us from the collinear divergence of dipole emission is unphysical, and as such, any conclusion obtained within our model must be strictly independent of its value. This means that we need to choose a value of r UV which is small enough to yield no effect on the shape of P n . At the same time, since the typical number of dipoles grows with r UV as a power, we can not pick a too small value, for the sake of saving computation time. In Fig. 11, we display how the multiplicity distribution gets altered if we vary the UV cutoff by a factor 2 , all other parameters being fixed:ᾱy = 4, r s /R = 0.025 and x 01 /R = 0.5. We find that the shape is unaffected by this choice. Simply, as could be expected, the number of dipoles is globally slightly lower for larger cutoffs. For the practical calculations, we have always set r UV = r s /100, which, we have tested, provides stable shapes of P n atᾱy = 4.