Hot spots and gluon field fluctuations as causes of eccentricity in small systems

We calculate eccentricities in high energy proton-nucleus collisions, by calculating correlation functions of the energy density field of the Glasma immediately after the collision event at proper time tau = 0. We separately consider the effects of color charge and geometrical hot spot fluctuations, analytically performing the averages over both in a dilute-dense limit. We show that geometric fluctuations of hot spots inside the proton are the dominant source of eccentricity whereas color charge fluctuations only give a negligible correction. The size and number of hot spots are the most important parameters characterizing the eccentricities.


I. INTRODUCTION
Collective azimuthal correlations, commonly parametrized in terms of harmonic "flow" coefficients v n of produced particles, were a crucial experimental signal in the discovery of the strongly interacting Quark-Gluon Plasma (QGP) in high-energy heavy ion collisions. More recently, very similar signals of collective correlations have also been discovered in high multiplicity events of various smaller collision systems, including proton-nucleus (pA) and proton-proton (pp) collisions at the Large Hadron Collider (LHC) [1][2][3][4][5][6][7][8][9][10][11], as well as proton/deuteron/helium-nucleus collisions at the Relativistic Heavy Ion Collider (RHIC) [12]. Since the lifetime of any such smaller collision system is significantly shorter than that of a heavy-ion collision, the emergence of such correlations was initially unexpected. This has led to an intense discussion [13] on whether the observed correlations in small collision system should be attributed to momentum-space correlations already present in the colliding projectiles [14][15][16][17][18][19][20][21], or whether they result from the final state response to the coordinate space geometry, either via hydrodynamical evolution [22][23][24] or by a simpler scattering mechanism [25][26][27][28][29][30][31].
Specifically, for the final state response to the coordinate space geometry, it was quickly realized that the subnucleonic degrees of freedom are crucial for understanding azimuthal correlations in pp and pA collisions [32][33][34][35][36][37][38][39]. This has lead to a series of new investigations into the transverse spatial distribution of subnucleonic degrees of freedom in protons and nuclei [34,[40][41][42][43][44]. Potential sources of subnucleon scale correlations in the energy density include both hot spot-like correlations, perhaps originating in the valence quark structure of the nucleon, and color charge fluctuations that ultimately result from the quantum mechanical randomness in the gluon radiation that generates the small-x degrees of freedom in a proton or nucleus. While a systematical description including both of these aspects is only beginning to take shape [45], both aspects of subnucleon scale correlations can already be included in specific models, and the purpose of this paper is to construct a simple and trans-parent model for the initial energy density in a protonnucleus collision, including some of the most important sources of fluctuations at a subnucleonic scale.
We will adopt the Color-Glass Condensate (CGC) [46] picture as a natural framework for discussing different sources of correlations and fluctuations of the initial energy density as a function of the transverse coordinate. We will include color charge fluctuations in a Gaussian model [47][48][49] parametrization in both the dilute probe and the dense target. The target nucleus will be treated as homogenous on average, with the Gaussian color charge fluctuations included to all orders in the color field. The inhomogenous probe proton, on the other hand, will be treated as dilute enough to linearize the energy density in its color charge, which also fluctuates as a Gaussian local charge density as in the MV model [47][48][49]. In addition to color charge fluctuations, we will consider the probe as consisting of N q hot spots of size r, with locations that fluctuate within the size of the probe R on an event-by-event basis. If these hot spots are thought of as originating from valence quarks, one would take N q = 3 for the proton, but we will remain more agnostic and consider N q as a free parameter, as in e.g. Ref. [38].
Notably, the degrees of freedom in our model are basically the same as in the IP-Glasma model [37,50,51], and to a first approximation our model could be considered as an analytically tractable simplification of IP-Glasma in the dilute-dense limit. However, two significant simplifications with respect to the IP-Glasma model need to be mentioned. First of all, unlike in IP-Glasma, we do not include any nucleon-level ("MC Glauber") structure on the nucleus side, but our target nucleus is always a homogenous (on average) and infinite sheet, corresponding to the large nuclear mass number (A → ∞) limit of the IP-Glasma model for the nucleus. The fluctuating nucleon positions inside the nucleus that are neglected here would naturally be important in order to describe the eccentricities of the created system. Secondly, while the Gaussian fluctuations of the color charges introduce some event-by-event fluctuations in the total energy produced, the IP-Glasma model also features additional "Q s fluctuations", which are not included in our model. In arXiv:2101.03791v1 [hep-ph] 11 Jan 2021 fact, the effect of such fluctuations is to make the color charge fluctuations non-Gaussian [19]. In addition to the eccentricities, these two sources of fluctuations are also important to describe the fluctuations in total multiplicity, i.e. different "centrality" classes. In fact, describing multiplicity fluctuations in proton-proton collisions has been a part of the motivation for introducing the "Q s fluctuations" [52,53].
Due to the very asymmetric treatment of the probe and target, our model might not be the best to fit experimental data, although we believe that it should be a good approximation for purely fluctuation-driven systems. Of course, in the future it could also extended by additional features such as nucleon position fluctuations in the target; however this is not the primary goal of this study. Instead, our purpose is to construct an analytically tractable model that enables us to cleanly study the effect of different sources of fluctuations, and of the values of parameters that describe them. The extreme dilute-dense limit enables us to perform all the averages over different types of fluctuations analytically. Since color charge fluctuations dominate on the side of the target, and hot spot positions fluctuations on the side of the probe, we can cleanly distinguish the effects of the two separately.
Our model depends on a set of phenomenological parameters that will be discussed explicitly below. These are the number of hot spots N q , the size of a hot spot r, the size of the proton R, and IR regulator for the long range Coulomb tails of the color field m, a UV cutoff C 0 and the saturation scales for the protonμ and the nucleus Q s . The same parameters are also present in the IP-Glasma model (where the lattice spacing a provides the UV cutoff). We will calculate initial energy densities at τ = 0, where they can be expressed analytically. We then obtain analytical expressions for the two-point energy density correlation function, which is the fundamental object characterizing the fluctuations and correlations of the energy density. From the energy density correlator we can then calculate eccentricities in a straightforward manner. The eccentricities are mostly driven by the long distance behavior of the energy density correlator. In particular this means that they are, contrary to the value of the energy density, quite insensitive to the UV cutoff. Since we use the two-point function of the energy density, we can only get "2-particle" eccentricities ε n {2}.
We will start this paper by briefly reviewing in Sec. II how the initial energy density in the Glasma is calculated from the Wilson lines representing the color fields of the colliding projectiles. We will then in Sec. III discuss how Wilson line correlators are calculated in our model for the nucleus, which is a the homogenous, infinite, fluctuating system of Gaussian color charges that is treated to all orders in the color fields. Further details on the Gaussian averaging procedure are given in Appendix A. We then introduce in Sec. IV our hot spot model for the proton, and analytically perform the geometric averages over the locations of the hot spots and over the color charges in the individual hot spots, which are treated as dilute objects. We then use this model to calculate the two-point correlation function of the energy density in Sec. V, demonstrating the sensitivity (or insensitivity) of these quantities to the parameters in our model. We then continue in Sec. VI to obtain values of the n = 2, 3, 4 harmonics, i.e. the eccentricities ε n {2}. Further checks of the dependence on the parameters are made in Appendix B. For comparison, some results are worked out in Appendix C in a model similar in spirit to that introduced in Ref. [54], where the energy density is a superposition of purely pointlike hot spots, which are only correlated through the constraint that their center of mass lies at the origin. We then conclude in Sec. VII, pointing out future directions and discussing how our model could be further constrained in the future.

II. ENERGY DENSITY IN GLASMA
In the CGC effective theory a high energy nucleon or nucleus is described as a cloud of small-x gluons radiated by large-x partons, represented as color charges. The gluon field is taken to be so dense that it can be described by a classical color field. The color field of a lone nucleus can then be expressed analytically as a function of its color charges. When the result for the color field is gauge rotated to the light cone gauge, one finds that the field is purely transverse and the transverse components are gauge rotations of the vacuum, sometimes referred to as "transverse pure gauge fields". The expressions for the color fields of two nuclei far away from each other, in their respective light cone gauges, are [48,49] Here U x and V x are light-like Wilson lines, the relevant degrees of freedom in CGC, that describe the eikonal interaction of a color charge moving through a color field. In a heavy ion collision one has two such color fields passing through and interacting with each other forming a non-equilibrium Glasma state [55], which continues to evolve into a Quark-Gluon Plasma. Immediately after the collision the interaction between the gluon clouds of the colliding CGC sheets leads to the creation of a longitudinal component of the gauge field. The gauge potentials can be expressed, in Fock-Schwinger (A τ = 0) gauge and at proper time τ = 0 + , as [56,57] A i where η is the space-time rapidity. Based on the results in Eq. (2) one can compute the energy density at τ = 0 + [55,58,59] and straightforwardly get the energy density two-point function Here we have expanded the gluon fields in the Lie algebra The brackets . in Eqs. (3) and (4) refer to an averaging over the color charges, or distributions of Wilson lines, in the two colliding nuclei separately. Since the color fields of the nuclei are built up slowly a long time before the collision, the two nuclei are not correlated. Thus the expectation values for the colliding nuclei factorize from each other and can be computed separately. In this work, we take the α-part to describe a dense nucleus and the β-part to describe a dilute proton. We will also separate the energy density two-point function into distinct connected and disconnected parts. We introduce the following decomposition of the four α and β correlators: where disconnected contributions α i,a x α k,c x α i ,a y α k ,c y refer to the parts appearing in expectation value of the energy density, while the connected part α i,a x α k,c x α i ,a y α k ,c y C of a correlator of α's or β's refers to the full correlator minus the disconnected part. The disconnected part of the energy density two-point function is just the product of energy density expectation values, given by the product of the disconnected parts of the α-and β-correlators. We refer to the contribution α i,a x α k,c y C from the disconnected part of the nucleus (α) correlator and the connected part of the proton (β) correlator as the "proton fluctuation" part of the energy density correlator. Analogously the connected-nucleus, disconnected-proton contribution α i,a x α k,c y is referred to as the "nucleus fluctuation" part of the energy density two-point function.
We expect that in our approximation the fully connected contribution (on both proton and nucleus sides) only gives a small contribution as it is sensitive to the fluctuations of both the proton and the nucleus, and we will neglect it in our calculation of the energy density two-point function in Eq. (4).
We will now go on to compute the 4-α correlator for the nucleus in the Sec. III and the 4-β correlator for the proton in the Sec. IV. The nucleus part will be computed as a full nonlinear Gaussian and the proton will be expanded to the first order in the proton saturation scale in the dilute limit.

III. NUCLEUS: FUNDAMENTAL REPRESENTATION 8-POINT CORRELATORS
Now we will compute the 4-point function of the nucleus gluon fields α as a nonlinear Gaussian expectation value, where we use the term "Gaussian" in the sense of Gaussian contractions of color charge densities. The calculation yields a result in terms of a generic two-point function of Wilson lines. We will then for simplicity adopt for the nucleus the GBW [60] parametrization for this two-point function. The GBW form is Gaussian in another sense, namely that it assumes a Gaussian functional form for the Wilson line dipole expectation value as a function of the transverse coordinate separation The GBW parametrization yields slightly simpler final expressions, since the crossed partial derivatives ∂ 1 x ∂ 2 x of the pure second power (x − y) 2 in the exponent vanish 1 . The algorithm for calculating Wilson line expectation values with Gaussian color charges is discussed in more detail in Appendix A. Here the nuclear saturation scale Q s is the only parameter needed to characterize our infinite, homogenous target nucleus.
We start by expanding the α's using their definition in Eq. (5), yielding Next we introduce new transverse coordinates to be able to pull out the derivatives and we write the traces in index notation, such that 1 In other words the linearly polarized gluon distribution vanishes [58].
This 8-point function of Wilson lines can be computed, using the algorithm used in [61], as Here the 24×24 matrix M 24×24 depends on the two-point correlator of Wilson lines at the coordinates x i , y i and color factors, and can be obtained as discussed in Appendix A. We will not write the rather lengthy expression here. The color structures δ a1a2 δ a3a4 δ a5a6 δ a7a8 etc. on the left are the elements of the chosen basis for different singlet operators that can be formed from the 8 Wilson lines, and the structure t a a2a1 t c a4a3 t a a6a5 t c a8a7 on the right corresponds to the particular correlator that we need to calculate here.
Instead of attempting to explicitly exponentiate the full 24×24 matrix, it is useful to first take the derivatives. We use the following identity for the derivative of a matrix exponential where M ≡ M n×n . As one can see, the derivatives do not enter the matrix exponentials themselves. Thus after taking the derivatives, we can take the coordinate limits x i → x, y i → y in the matrix exponentials, rendering their evaluation easier. The result is not written here for it is a long equation and is not a new result. It was computed with a slightly different formulation of the same algorithm in Ref. [59] and it corresponds to the result we get using this method.

IV. PROTON: HOT SPOT MODEL
We will treat the proton as a collection of Gaussian hot spots. The hot spots consist of color charges with Gaussian fluctuations that are local (uncorrelated between different points), i.e. given by the MV model [47][48][49]. The hot spot approach is inspired by the picture of a proton having three large-x valence quarks and has been used in the IP-Glasma model [40,41]. We also let these hot spots to be Gaussianly distributed in the transverse coordinate inside the proton with respect to its center of mass.
We now have two averages to calculate in order to obtain the energy density correlator. We define a double average of an operator (12) to include both the average over the color sources in the MV-model, and the averaging over the hot spot coordinates. Since we eventually want to calculate eccentricities with respect to the center of proton, which fluctuates with the hot spot positions, it is essential to explicitly fix the center of mass of the hot spot system 2 to a known coordinate B. The distribution of the hot spot locations is taken to be a Gaussian: where the parameter R has an interpretation as the proton radius. The prefactor in Eq. (12) is chosen so that the expectation value is normalized: 1 = 1.
The distribution of pointlike color charges within a hot spot is taken to have a Gaussian distribution where r is the radius of the hot spot and µ 0 is a coefficient characterizing the color charge density. The hot spots enter the calculation through the two-point function of color charge densities as follows where N q is the number of hot spots in the proton. Correlators of more than two color charge densities are calculated by taking the distribution to be Gaussian, and can thus be expressed in terms of the two point function (15). Now we can compute the proton contribution to the energy density two-point function.
which we will do in the limit of small charge density, i.e. to lowest order in the parameter µ 0 . We start by writing the color fields β in terms of the Wilson lines We then expand to lowest order in the color sources yielding which relates the color field in the Wilson line to the color charge density. We have regularized the infrared behavior of the Green's function with a mass m, which should be thought of as a confinement scale regulator , can, assuming x = y, be written as where K 1 is the modified Bessel function of the second kind. This expression contains an ultraviolet, shortdistance singularity at points with x = y, which arises from the fact that our color charges are treated as a pointlike objects. We regularize this by a short-distance cutoff, making the substitution where Θ is the Heaviside step function and C 0 is a short distance cutoff. Now the color charge and hot spot averages can be evaluated fully, giving the correlator as This is explicitly factorized into a part describing the coordinates of the color charges a, b, and the Green's function part describing the color field generated by these charges. In Eq. (22) the averages over the color charges have naturally split into two kinds of contributions. The first one, proportional to the number of hot spots N q , results from taking all color charges from the same hot spot, resulting in the function where the second term in the exponent makes it manifest that the color charge coordinates a, b are within a distance ∼ r from each other. The second contribution results from taking color charges from two separate hot spots, and is proportional to N q (N q − 1), the number of pairs of distinct hot spots. It is given by the function where i, j ∈ {1, . . . , N q } and i = j. Here it is also clear that the coordinates of the color charges a, b are typically separated by a distance ∼ R of the order of the size of the proton. Let us finally summarize here the phenomenological parameters characterizing our hot spot model. They are the proton radius parameter R introduced in Eq. (13), the proton color charge density µ 0 and hot spot size r introduced in Eq. (14), the number of hot spots N q introduced in Eqs. (12), (15), the "gluon mass" infrared regulator for the Coulomb tails of the color field m in Eq. (19), and the short distance cutoff C 0 in Eq. (21). We will discuss the dependence of the energy density correlator and the eccentricities on these parameters in the following sections.

V. RESULTS: ENERGY DENSITY AND ITS CORRELATOR
Let us start by evaluating the energy density expectation value (3) in the framework of a full nonlinear Gaussian nucleus and hot spot model proton. Let us first compute the part with the nucleus side average. We start by using the definition of the α's and writing the resulting expression in index notation. We will also introduce new transverse coordinates to enable us to pull the derivatives out of the correlator. Doing this we get Again, using the algorithm presented in [61] and discussed in more detail in Appendix A, we can express this as This reduces to Now we can take the derivatives of the matrix exponential and then take the coordinate limits in the similar fashion as in section III. Doing this, and using the GBW model as an input, one finds Plugging this back to the expression of the energy density, Eq. (3), and doing the color and transverse index algebra, we get where N is the number of colors. Now we will calculate the proton side contribution in the hot spot model in the dilute limit. We will start by using the definition the β's and expanding the Wilson lines to the lowest order in sources yielding Taking the CGC and hot spot averages, we get where we defined, analogously to Eqs. (23) and (24), with i ∈ {1, . . . , N q }. The function F 1 can be interpreted as the average density of proton color charges in the transverse plane.
Now that we have obtained the expression (31) for the one-point function of the energy density, we move to calculating the two point function, i.e. the energy density correlator. We want to separate the contributions due to the proton and nucleus fluctuations. To do this, we consider separately the disconnected and connected part of the proton and nucleus contributions as defined in Eq. (6). The disconnected-disconnected contribution of the energy density two-point function is just the product of two energy densities in the two different transverse positions Note that this expression is radially symmetric in both transverse coordinates and can not have an explicit correlation between the two coordinates.
Next we want to compute the nucleus disconnected and proton connected part, i.e. the "proton fluctuation" contribution. Subtracting the disconnected-disconnected contribution does not lead to a particularly simple expression (see more detailed discussion in Appendix B), so we write this part here as with the same F 2 and F 3 as in equations (23) and (24). Lastly, we need the nucleus connected and proton disconnected or the "nucleus fluctuation" contribution, which takes the form As discussed earlier, we will neglect the connectedconnected contribution here. Thus we now have the results for the energy density correlator that we will need. Let us now evaluate numerically the different parts of the energy density two-point function. We will plot both the energy density correlator values themselves, and below these the relative contributions of the different parts that contribute to the total two-point function (in our approximation where we drop the connected-connected part). We will vary parameters of our model to produce an error band for the energy density correlators, with the relative contributions only plotted for the central values in this variation for clarity. We have normalized the twopoint functions in the top plots with the nucleus saturation scale Q s and the scaled proton saturation scale parameterμ 2 defined as In the rest of the paper, we use the following default parameters for plots unless stated otherwise. We always have three colors N c ≡ N = 3. The nucleus saturation scale is set to be Q s = 2 GeV, the UV cutoff is C 0 = 0.05 GeV −1 , the IR regulating mass is m = 0.22 GeV, the proton radius parameter is R = √ 3.3 GeV −1 and the hot spot radius is r = √ 0.7 GeV −1 . For this work, the values for R and r were taken from [42]. The mass was chosen to be of the order of the QCD scale and C 0 was taken to be of a scale smaller than the hot spot radius to allow for gluon fields originating from the same hot spot to overlap but still cut off the divergent behavior of the Green's functions.
We will plot the energy density correlator in two different coordinate configurations. Firstly, we set the two coordinates in the two-point function to be equal (x = y) and plot the parts of the two-point function as a function of the distance from the center of the proton (B) divided by the proton radius R. This measures the local energy density fluctuations as a function of distance. These plots are shown in figure 1. Secondly we have plots of the energy density correlation, where we take a straight line through the center of the proton and let the two coordinates move, at the same rate, in opposite directions along the line. In the plots the two-point function parts are plotted as functions of the distance of these two coordinates divided by the proton radius. These plots are shown in figure 2.
We add error bands to the plots by varying either the UV cutoff (C 0 ) or the IR regulator (m) by ±50%. The UV cutoff dependence seems to be much more prominent than the dependence on the IR regulator. This is not surprising, it has been long known that the energy density at exactly τ = 0 is logarithmically UV divergent. This divergence would, however, disappear with the evolution to larger τ [58,62], which we do not attempt to do here 3 . The UV-cutoff dependence is particularly clearly visible near the center of the proton, where the color charge density is largest and there is large amount of short range overlap for the gluon fields generated by these color charges. The mass, on the other hand, has more of an influence on the long range behavior on the gluon fields. The long range mass dependence is shown in figure 3, which indeed shows that the exact value of the mass influences where the proton fluctuations are sizable or larger than the contribution from the fully disconnected part. We further note that the contribution from proton fluctuations to the energy density correlation, Eq. (34), is actually a sum of two contributions, a short range one originating from a single hot spot (F 2 ) and another one that is sensitive to two hot spots (F 3 ). The latter is preferentially long range, because the hot spots are typically separated by a distance ∼ R. The interplay between these two contributions leads to an interesting structure of the correlation function with a dip and secondary maximum seen in Fig. 3. These contributions are shown separately in Fig. 8 in Appendix B.
The dependence on the number of hot spots (N q ) can be seen in Figs. 1 and 2. When we only have one hot spot, in our model, it has to be in the center of the proton, and thus there are only color charge fluctuations for the proton. This results in a small contribution for both the proton and the nucleus fluctuations. With a small, larger than one, number of hot spots we start to see that, in some regions, the proton fluctuations can even become larger than the fully disconnected contribution. Especially the largest fluctuations at the edge of the proton seem to be largely due to proton shape fluctuations. In the plots where we increase the separation of the two coordinates, we can see that at small separation, the proton fluctuations are sizable in comparison to the fully disconnected contribution and at large separation they can become even larger than the disconnected contribution. However, this behavior largely disappears with a large number of hot spots. The nucleus fluctuations are Overall, we see that the hot spot fluctuations are by far the dominant mechanism for fluctuations and correlations of the energy density. While the exact magnitude of their contribution depends on the values of the parameters, the numbers used here correspond quite well to the ones currently used in phenomenology. The value of the energy density depends strongly on the UV-cutoff, as expected. The dependence on the IR-regulator m is less pronounced, but it is important in the long distance tail of the distribution. Consequently, it turns out to be very significant for the eccentricities, as we will see in the following Section.

VI. RESULTS: ECCENTRICITIES
Now let us move to computing the azimuthal eccentricities arising from the proton shape and gluon fluctuations in our model.
For a single configuration, the deformation of the energy density field ε(x) from azimuthal symmetry is quantified by the dimensionless ratios known as the eccentric- ities, defined as [54] where B is the center of the event defined as Event-by-event Monte-Carlo studies in IP-Glasma determine the center of mass B from the energy density as in Eq. (38). Here, in contrast, we will neglect the effect of the color charge fluctuations on the location of the center of the system, and take the center of the collision system as the center of mass of the hot spots, as defined in Eq. (12). This is justified by the assumption that the hot spot location fluctuations are much larger than the gluon field fluctuations. The fluctuation part of the energy density is defined for individual events as The average of the fluctuations δε(x) vanishes. Thus to actually have sensitivity to the fluctuations in a fluctuation driven system, we need consider a different quantity that is somehow quadratic in the energy density. One could calculate the eccentricity from the mean square of the ratio (37) as Calculating the expectation value of the ratio is, however, not desirable for several reasons (see also the related discussion in Ref. [21]). Firstly it is not possible to evaluate (40) based on just the two point function of the energy density that we have available in our model. Nevertheless we consider the energy density correlator as the more fundamental quantity that should contain the relevant physics. We could only evaluate (40) in the regime of small fluctuations, but this might not be a valid approximation for us due to the large fluctuations arising from the hot spot model. More generally in the CGC formalism, or in quantum mechanics in general, it is natural to calculate expectation values of operators that are positive integer powers of physical observables, such as the energy density. Also in experimental determinations of flow coefficients v n , it is not immediately obvious to us whether one is closer to calculating ratios of expectation values or expectation values of ratios; this depends on how the event-by-event varying total multiplicity is treated in the analysis. In light of this discussion we will here define the eccentricity through the more natural quantity, namely the one obtained from the ratio of the averages of azimuthal harmonics of the energy density two point function. The definition of the eccentricity that we use is Due to the infinitely large homogeneous nucleus in our model, this quantity is actually independent of B, and we can set B = 0 without loss of generality. We note that our definition in Eq. (41) is equivalent to the expectation value of the ratio (40) when the fluctuations of energy density (39) are small. In this case one can expand the eccentricities to the first order in the fluctuations, which leads to the approximation used e.g. in Ref. [54] where the difference with respect to Eq. (41) is that the two-point function in the denominator gets replaced by the product of one-point functions. This approximation has been used in [65] in the case of A-A collisions when considering only the color fluctuations of the nuclei. Now that we incorporate the shape fluctuations to the proton through the hot spot model resulting in larger fluctuations, we can not a priori be sure that the use of this small fluctuation approximation is valid in our case. However, in practice the approximation (42) is quite close to our full result, as demonstrated in Appendix B.
Our main results for the eccentricities ε 2 {2}, ε 3 {2} and ε 4 {2} are shown in Figs. 4, 5 and 6. We compute the contribution to the eccentricities separately for proton and nucleus side fluctuations, noting that the sum of the proton and nucleus side fluctuation eccentricities does not equal to the total eccentricity, since the contributions are added quadratically (c.f. Eq. (41)). Figure 4 shows the eccentricities for n = 2, 3, 4 for the default values of our parameters with three hot spots (N q = 3). We see that the proton fluctuation contribution completely dominates the eccentricities. The error bars, resulting from a variation of the parameter m by ±50%, show that the result is very strongly dependent on this infrared cutoff, especially for the higher harmonics. It is important to to emphasize that, in contrast to the impression one would get from the two-point energy density plots, the dependence on the UV and IR regulators has changed. The UV regulator mainly affects the normalization of the energy density correlator, but has very little influence on its coordinate dependence. Thus it has very little effect on the eccentricities, as quantified explicitly in Fig. 7 in Appendix B. Conversely the IR regulator affects the long range physics to which the eccentricities are sensitive due to the large ∼ |x − B| n weight of the long range tails, resulting in the sizeable uncertainties of the eccentricities seen in Fig. 4.
We then focus on the dependence of the eccentricities on the number of hot spots, which is plotted in Fig. 5. One can see that the overall trend for the proton fluctuation part of the eccentricity is such that it decreases as the number of hot spots increases. This can be explained by the reduction of event-by-event fluctuation due to a large amount of independent hot spots, which makes the events smoother and more azimuthally symmetric. We also see that the nucleus fluctuations are not as sensitive to the number of hot spots, as expected. The proton fluctuations remain the dominant source of eccentricity up to very high numbers of hot spots, and the nucleus fluctuation contribution is negligible. An important exception to this systematic trend is ε 3 at N q = 2, where the hot spots are forced to be exactly back-to-back, and the contribution to ε 3 is solely due to color charge fluctuations.
Next we study the effect of the size r of an individual hot spot. Figure 6 shows the second order (n = 2) eccentricity as a function of the hot spot size r divided by the proton size parameter R (obtained keeping R constant and varying r). Firstly we notice that the dependence on the mass is still rather large in comparison to the ex-act value of the hot spot radius. The hot spot radius does have a sizable effect on the proton fluctuation part, with larger hot spots making the system smoother and consequently less eccentric, as expected.
Overall, we find that the values of the eccentricities can be substantial, with ε 2 {2} reaching above 0.5 for a small number of hot spots, with the exact value very strongly depending on m. In order to set these values in perspective, we have also constructed a model of an energy density constituted by N q exactly pointlike hot spots, correlated only by the requirement that their center of mass is at a fixed coordinate. This model, inspired by but not equivalent to the work in Refs. [54,65], would correspond to the limit r → 0, m → ∞ of our hot spot model, where energy density consists of a superposition of delta function like peaks. In this limit, the eccentricities can in fact be calculated analytically, and the results are provided in Appendix C. Figure 5 also shows a comparison of the eccentricities to this pointlike energy density model. One sees that the eccentricities of the hot spot model are still a factor of ∼ 2 smaller. We therefore conclude that the smoothing provided by a finite hot spot size r and the Coulomb tails of the color fields regulated by m, still has a significant impact on the eccentricities, and a realistic color field calculation cannot be accurately described by an energy density consisting of just a finite number of delta function peaks.

VII. CONCLUSIONS
We have constructed an analytical model for calculating the initial energy density and its fluctuations in the initial stages of a high energy hadronic collision. Our model is formulated in an extremely asymmetric dilutedense limit, which enables an analytical calculation of the energy density correlator. In this limit, we can very explicitly see the influence of different physical aspects of the model (a philosophy very strongly advocated recently e.g. in Ref. [66]). The most important aspect that we have focused on here is the role of hot spot fluctuations in the dilute proton projectile, which have recently been understood to be crucial for understanding the geometry of the initial stages of high energy collisions [39]. We have also quantified the extremely important role for the eccentricities played by the infrared regulator parameter m. In contrast, the UV divergence present in the glasma fields in the MV model at τ = 0 has very little influence on the eccentricities.
In our model, with parameters taken to have what we believe as reasonable values for a central proton-nucleus collision, the proton hot spot fluctuations are much more important for the eccentricities than the color field fluctuations in the nucleus. In addition to the IR regulator for the Coulomb tails m, their effect strongly depends on the hot spot size r, and the number of hot spots N q , which we treat as a free parameter. In order to enable a more realistic comparison to experimental data in the fu-ture, it would be important to include geometric (nucleon position and internal nucleon structure) fluctuations also in the target nucleus. We have not done so here, to be able to more cleanly demonstrate the effect of the separate sources of fluctuations. Another future development is the use diffractive ep scattering data to constrain the parameters of our model. For example, we believe there is a significant degeneracy in the numerical values of the hot spot size r and the IR cutoff m. Our model for the proton color fields is simultaneously versatile enough to describe the relevant physics, but simple enough to be analytically solvable. We believe that this combination will make it useful for a systematic study of such effects, which we plan to pursue in a future publication.
Another future avenue of study wold be to combine this initial condition with an analytical approach for the τ dependence, e.g. through some suitable resummation of an approach like the one of Refs. [63,64]. It would also be interesting to extend our calculation to 4th order correlators of the energy density. This would be possible in principle in our framework, but especially on the nucleus side the calculations would become rather complicated as the required 16-point fundamental Wilson line correlators have never been calculated so far.
Acknowledgements: We thank A. Soto Ontoso, B. Schenke, H. Mäntysaari and P. Guerrero Rodríguez for insightful discussions. We gratefully acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC-TR 211 "Strong-interaction matter under extreme conditions" Project number 315477589 (S.S.) and the Academy of Finland, project No. 321840 (TL). This work is supported by the European Research Council, grant ERC-2015-CoG-681707. The content of this article does not reflect the official opinion of the European Union and responsibility for the information and views expressed therein lies entirely with the authors.

Appendix A: Computation of the nucleus side Wilson line correlators
We wish to evaluate the nucleus side Wilson line correlators in a Gaussian model for an infinite homogenous nucleus. The Wilson lines can be written as and are defined through the series expansion of the exponential. In this Appendix we follow the method presented in [61] where it was used to compute 4-point Wilson line correlators.
As a first illustration of the method let us discuss the simpler case of computing the 2-α correlator needed for the energy density and the nucleus disconnected part. First we have to choose a basis for the computation. The exact choice depends on the correlator we are computing. In our case we can choose the basis to be built from basic building blocks consisting of two Wilson lines with open fundamental indices Now the actual correlator we wish to evaluate is The main idea in the computation is to expand the Wilson line exponentials into series and then taking all the possible contractions of two color charge densities ρ, as permitted by the Gaussian model for the CGC weighting functional. The main ingredient for the calculation, the 2-point function of ρ, reads where µ 2 (x + ) is the density of color charges. One also has to take into account the path-ordering of the Wilson lines, which results in all but two types of contractions to give a zero. The tadpole type contractions factorize out, and can be recombined to the final result later. On the other hand, the contractions that connect two different Wilson lines result in state transitions and are the part that is nontrivial to compute. Organizing the calculation in such a way that we start to contract the Wilson lines starting from LC time z + = ∞ and progressing towards z + = −∞, the state transitions do not care about how the Wilson lines have been contracted in the distant past. In fundamental representation contractions, which we are right now dealing with, the color algebraic structure of the states changes according to the well known Fierz identity t a ij t a kl = Considering this, we can see that the states that mix with our correlator are the ones where we cut open the z + = ∞ side of the Wilson line "dipoles" and connect them in every possible way while still connecting only one Wilson line to one daggered Wilson line. Thus with 4 Wilson lines we have two possible basis states which can be chosen to be Let us now express a general correlator in this basis as a vector, where standard basis vectors represent the different basis states in such a way that W i → e i . In this way, the general correlator can be written as Next one needs to compute the transition matrix that describes how the basis states evolve as one contracts two ρs, each from a different Wilson line, in every possible way. This matrix gains more and more powers as we contract more and more ρs. These powers of the matrix exponentiate and one gets an expression where the last nontrivial task is to compute the matrix exponential. The expression still requires initial conditions for the basis states, which we again express as a vector. This means that we need to evaluate the basis states when all the ρs have been contracted i.e. when the Wilson lines go to identity. In our case, they tell how the fundamental indices of the Wilson lines are connected in the end. For example the initial condition vector for our chosen basis states would be Now the full expression for the 4-point correlator we need for the nucleus can be written as The transitions 1 0 The correlator , (A9) where the matrix M 2×2 is proportional to the density that has been integrated over the longitudinal coordinate The expression (A9) for the correlator can be broken into 3 parts: the initial conditions and the transition matrix, which are universal for any correlator that can be expressed in our chosen basis, and the correlator vector, which depends on the actual correlator we wish to compute. Now that the expression we are dealing with has fundamental representation SU (N ) group generators contracting the correlators, we can just contract both sides of this equation with the t a k2k1 t a k4k3 . Doing this we get which is exactly the correlator we needed for the 2-α correlator For details for finding the transition matrix for a given basis see [61], where it has been done for the 4-point Wilson line correlator. We present an algorithm for finding the matrix for the specific basis we have chosen here. Firstly let us give a name for the basis states for easier handling. Let us denote them with an S with 2-tuples as arguments, which tell what fundamental index is associated with which Wilson line and which Wilson lines are connected to each other. Each of these "open dipoles" are separated by a semicolon. For example we can write the first basis state as The other basis state is just a permutation of the 2-tuples in such a way that we pair the Wilson lines with the daggered Wilson lines in a different way.
To find the transition matrix, we compute how the basis states evolve in one step of evolution i.e. when we have one fundamental color matrix contraction done in every possible way, omitting tadpoles for now. When organizing the computation in this way, it suffices to compute the evolution of one basis state as the other differs from it just by a permutation of 2-tuples. The evolution for the first basis state is In the MV model, one would have where the Green's function is the same as the one defined in (19). Now by identifying the basis states and making the replacement W i → e i , one has found a column vector of the transition matrix. If one does the evolution in equation (A14) for the i:th basis state, the resulting column vector is the i:th row of the transition matrix. Doing this for all the basis states, one finds the transition matrix, which is still missing the tadpole contributions. The tadpole contributions always combine with the ladder-type contributions in the same way forming the quantities which can be computed in the MV model using the definition of the L-function. The final transition matrix can then in practice be found with the replacement As an example, the resulting transition matrix M 2×2 , with this choice of basis, is . (A18) Here we have used the short hand notation In this work we use the GBW parametrization to describe the nucleus gluon fields. In Gaussian models the Wilson line correlators can be expressed in terms of the Wilson line dipole. One can thus match the Γ-function appearing in the dipole to the GBW dipole and use this Γ to obtain the higher point correlators. By computing the Wilson line dipole using the method we have used for the other correlators in this work, one can find that it can be written as On the other hand the GBW dipole can be written as [60] D GBW By matching these, one finds that the required transformation for expressing our results in the GBW model is The same calculation can be done for the fundamental 8-point Wilson line correlator that is required for the nucleus connected part of the energy 2-point function. In this case we have 24 states that mix with our correlator. Again, the correlators can be built from the basic building blocks A2 and can be found by pairing the Wilson lines with the daggered Wilson lines in every possible way. The correlator we wish to compute is (A23) The exact choice of basis does not matter so we just set the state corresponding to the correlator to be the first basis state. Again, we write the basis states in a more handleable manner as The evolution of other basis states can be found by permuting the 2-tuples to match them. After this, we can again do the replacement W i → e i and use the result as the i:th column of the transition matrix M 24×24 . Again, the tadpole contribution can be added to the matrix by the replacement in Eq. (A17).
Now we have all we need to compute the 4-α correlator we needed for the nucleus connected contribution As the evaluation of the matrix exponential is the most cumbersome part of computing this correlator, one can use the identity to simplify the calculation. One can take the derivatives using the identity and then take the coordinate limits in the exponents, which causes some of the Γ-functions to vanish.
The expression we obtain for the color contracted nucleus correlator is We have assumed that the model used has the property but did not make further assumptions of the model used for the infinite, homogenous nucleus. The correlator in Eq. (A28) is the one needed for the contraction of the proton disconnected part with the nucleus connected part of the energy density 2-point function. The superscripts on the Γ-functions indicate derivatives. i is with respect to x 2 , j is with respect to x 4 , k is with respect to x 6 and l is with respect to x 8 . So for example After taking all the derivatives, one has to take the limits x 1,2,3,4 → x and x 5,6,7,8 → y. We used the FeynCalc package to handle the required color algebra for this correlator [67][68][69].  Our model required two different regularization parameters, one for the UV and one for the IR divergences. When presenting results for the eccentricities, we have used mass variation to obtain confidence bands for our results, as that was the major source of uncertainty among these two parameters. In figure 7 we present some example plots that show that the UV cutoff (C 0 ) dependence of the eccentricities is indeed much smaller than their dependence on the IR mass regulator (m).  In figure 8 we have plotted the different parts of the proton connected contribution to the 2-point function of the energy density to explain the emergence of the second peak in the right-hand side plot in figure 3. Different contributions to the proton fluctuations can be separated into the following contributions F2 with disconnected Green's functions F2 with connected Green's functions F3 with disconnected Green's functions F3 with connected Green's functions (B1) where the disconnected contribution is subtracted proportionally from the terms with the disconnected Green's functions as it seems like the natural way to distribute it between the two contributions. In this way all of the disconnected contribution is subtracted from the F 2 term when N q = 1 i.e. when the F 3 term vanishes. Then an increasing fraction of the disconnected contribution is subtracted from the F 3 term as it grows with a growing N q . By inspecting Fig. 8 one observes that the first peak in the proton fluctuations is dominated by the contributions from a single hot spot (F 2 ), while the second peak emerges when the contribution from two hot spots (F 3 ) becomes dominant. In figure 9 we show how our definition of eccentricity (41) differs from the approximation in equation (42). In the parameter space we have been exploring in this paper, the results of these two definitions do not differ much. However, as one moves to the parts of the parameter space, where fluctuations are large, for example when the hot spots are very localized, one starts to see that the the result obtained from the two different definitions start to differ more. In the extreme case of energy density being a collection of Dirac delta-like hot spots as in Appendix C, the approximation (42) is not bounded by one unlike the definition (41) that we use throughout this paper.

Appendix C: Pointlike energy density model
We consider the limiting behaviour for our model by taking the energy density to be a collection of Dirac delta hot spots of extremely localized gluon fields. This should be the limit of our model if we formally took the limits of zero hot spot size (r → 0) and infinite mass (m → ∞). The expression for the energy density is now where, as before, b i denotes the position of the hot spots and ε 0 is a dimensionful constant whose exact value is of no importance as it cancels out in the end. Now having the expression for the energy density, we can compute its one and two point functions using the hot spot averaging procedure yielding which is valid for N q ≥ 2, and which is valid for N q ≥ 3.
Based on the above results of the one and two-point correlation functions of the energy density, we now we go on and compute eccentricities ε n {2} defined in Eq. (41) in this model. For the numerator we need to compute the weighted integral of the two-point function: d 2 xd 2 y|x − B| n |y − B| n cos(nθ x−B − nθ y−B ) ε(x)ε(y) = ε 2 0 R 2n (−2) n (N q − 1)N 1−n q nΓ(n) + ε 2 0 R 2n N q 2 − 2 N q n nΓ(n), (C4) where here and in subsequent equations Γ(n) is the usual gamma function. We further need the denominator integral for the two-point function, and to be able to compare the definitions of eccentricity, the integral of the energy density one-point function squared: for N q ≥ 3 and we expect these values to give the absolute upper limit for the eccentricities in our hot spot model. By looking at the results in Fig. 5 one can notice that in the point like energy density model, the n = 3 eccentricity is smaller than the n = 4 eccentricity. This is due to the fact that we have fixed the center of the proton and thus there is a preference for the hot spots to be "back-to-back". This effect shows up as the first term in (C4) being proportional to (−1) n . Now because of this, the relative angle between x and y, the two points we measure the energy density at, is going to preferably be close to π. Thus the cos(nθ x − nθ y ) is preferably going to be cos(nπ), which gives a positive sign for the n even eccentricities and a negative sign for the n odd eccentricities as seen in (C4). The other term is due to the gluon fields being taken from the same pointlike source, for which the cosine always gives a 1. A hint of this same behavior can be observed in our actual model when considering the upper error bands of our eccentricities with a small number of hot spots. This makes sense as the upper bound corresponds to a larger mass which results in a more localized gluon field originating from a hot spot.