Symmetric cumulants as a probe of the proton substructure at LHC energies

We present a systematic study of the normalized symmetric cumulants, NSC(n,m), at the eccentricity level in proton-proton interactions at $\sqrt s\!=\! 13$ TeV within a wounded hot spot approach. We focus our attention on the influence of spatial correlations between the proton constituents, in our case gluonic hot spots, on this observable. We notice that the presence of short-range repulsive correlations between the hot spots systematically decreases the values of NSC(2,3) and NSC(2,4) in mid-to-ultra central collisions while increases them in peripheral interactions. In the case of NSC(2,3) we find that, as suggested by data, an anti-correlation of $\varepsilon_2$ and $\varepsilon_3$ in ultra-central collisions, i.e. NSC(2,3)$<0$, is possible within the correlated scenario while it never occurs without correlations. We attribute this fact to the decisive role of correlations on enlarging the probability of interaction topologies that reduce the value of NSC(2,3) and, eventually, make it negative. Further, we explore the dependence of our conclusions on the values of the hot spot radius and the repulsive core distance. Our results add evidence to the idea that considering spatial correlations between the subnucleonic degrees of freedom of the proton may have a strong impact on the initial state properties of proton-proton interactions [1].


I. INTRODUCTION
The vast amount of high-precision data collected in the recent years at the Large Hadron Collider (LHC) has allowed to go beyond the analysis of event-averaged observables and explore higher order moments of their probability distributions. In the context of heavy ion collisions, a very recent tool to study the properties of the Quark Gluon Plasma (QGP) are the correlations between different flow harmonics v n i.e. the symmetric cumulants defined as [2,3] SC(n, m) = v 2 or in their normalized version that eliminate the dependence on the absolute magnitude of v n(m) . The measurement of Eq. 1 would be zero by definition, if the fluctuations of v n and v m were totally uncorrelated, in the same way as the Pearson's correlation coefficient. Instead, a positive value of Eq. 1 implies that an event with v n > v n would be more likely to have v m > v m . In particular, by measuring SC (2,3) we can gain information about initial state fluctuations whereas SC (2,4) is mostly sensitive to the strongly interacting medium properties [4]. * albacete@ugr.es † petersen@fias.uni-frankfurt.de ‡ ontoso@fias.uni-frankfurt.de The experimental study of symmetric cumulants was pioneered by the ALICE Collaboration [5] in Pb+Pb collisions at √ s NN = 2.76 TeV. Recently, the CMS Collaboration has performed an experimental analysis of the symmetric cumulants as a function of the multiplicity in the three collision systems available at the LHC: p+p, p+Pb and Pb+Pb [6]. The experimental results suggest a similar pattern across systems. Concretely, SC (2,4) is always positive although its multiplicity dependence varies from p+p to Pb+Pb. Further, NSC (2,4) is clearly modified when varying the system size. On the contrary, the sign of SC(2,3) is strongly multiplicitydependent: at low multiplicities SC(2,3) is found to be positive. However, it turns out to be negative for very high multiplicities, N offline trk > 60 in both p+Pb and Pb+Pb and N offline trk ∼ 100 in p+p. Moreover, NSC (2,3) in the high-multiplicity regime is found to have not only the same sign in the three collision systems but the same quantitative value as well. It should be noted that in the p+p case although the trend of NSC (2,3) to smaller values when increasing the multiplicity is visible, the systematic uncertainties make it compatible with zero. Relating the negative sign of NSC(2,3) at very high multiplicities with the initial geometry of the proton is the main goal of this work.
Previous to the symmetric cumulants study the similarities between high-multiplicity proton-proton interactions with p+Pb and Pb+Pb were already observed in the individual flow harmonic coefficients [7][8][9] and the enhanced production of multi-strange hadrons measurements [10]. All together, the experimental data is constantly reigniting the debate on whether collective effects, precedently attributed to the formation of QGP droplets, are being observed in small collision systems such as p+p and p+Pb [11][12][13][14].
From a theoretical point of view, the well established paradigm of either Monte Carlo Glauber [4] or Color Glass Condensate [15][16][17] initial conditions followed by viscous hydrodynamic evolution has successfully described the data on SC(n,m) in the Pb+Pb case. Turning to smaller systems, the situation is less conclusive. On the one hand, there have been attempts to describe the values of SC(n,m) in p+Pb by computing them in terms of eccentricities, i.e. replacing v n(m) by ε n(m) in Eqs. 1-2, within wounded quark models. These studies lead to the correct negative sign of SC(2,3) at high multiplicities but the magnitude is off [18][19][20]. Also without hydrodynamic evolution, qualitative agreement with the p+Pb results has been achieved in an initial state model where the partons in the projectile coherently scatter off the color fields in the heavy nuclear target [21,22]. Up to today we are unaware of any theoretical prediction for the values of SC(n,m) in p+p interactions at LHC energies, although results for RHIC energies were presented in [18].
A key ingredient in the computation of the symmetric cumulants in any of these frameworks is the parametrization of the initial geometry of the collision. Especially the smaller systems exhibit a high degree of sensitivity to the description of the proton structure and its fluctuations. The importance of considering subnucleonic degrees of freedom when describing the elliptic flow in p+Pb within the IP-Glasma framework was realized in [23]. In parallel, considering the proton constituents to be subjected to spatial correlations has been shown to have a substantial impact on the values of ε 2 and ε 3 [1] beside other features of hadronic interactions such as the hollowness effect [24,25]. This work constitutes the natural extension of our previous studies on the initial state properties of proton-proton collisions in terms of eccentricities [1] by exploring not only their mean but their fluctuations. As we shall explain along the manuscript we rely, for simplicity, on a geometrical picture of the collision. Therefore we use Monte Carlo Glauber [26] simulations, where the proton is composed by, in principle, 3 gluonic hot spots. By relating the hot spots to the gluon clouds radiated by the valence quarks considering the proton to be formed by 3 constituents becomes natural. However, the possibility of having a different number, N hs , that may account for other elements such as the large-x sea quarks is explored in this work. The centrality selection is done in terms of the entropy deposition as a proxy of particle production. We find that the inclusion of short-range repulsive correlations has a critical impact on the sign of NSC (2,3) in ultra central collisions. The net effect of the presence of correlations is to reduce the value of NSC(2,3) with respect to the uncorrelated scenario in the more central collisions and even push it to negative values in the highest centrality bins. An intuitive interpretation of this result is given by characterizing the topology of the interaction in terms of the number of wounded hot spots and of collisions between them. In the case of NSC (2,4), the results are qualitatively the same as for NSC (2,3) although it always remains positive within the regions of the parameter space explored in this work. In order to disentangle the possible phenomena that could contribute to the negative sign of NSC(2,3) we compute it for different values of the parameter space. We conclude that, as expected, not only the presence of correlations is important but also the interplay between the different scales of the problem, that is, the number of gluonic hot spots, their radius and the value of the repulsive distance. The radius of the hot spot, R hs , has been studied, in terms of the correlation length of the gluon field strengths inside hadrons, via lattice QCD calculations [27] and within perturbative [28] and non-perturbative [29] frameworks. On the contrary the value of the repulsive distance, r c , apart from being different from zero [25], is essentially unconstrained. So it is the number of gluonic hot spots. Thus, this study guided by the experimental data on NSC(n,m) helps to restrict the values of N hs , R hs and r c within our model.
The organization of the paper is as follow. We begin by reviewing the main ingredients of our model in Sec. II. Then, in Sec. III, the results for NSC(2,3) and NSC(2,4) as a function of centrality are presented. As the most interesting results occur on the 0-1% centrality bin, we focus on it in Sec. IV and study the role of the interaction topology together with a scan of the parameter space. Further, we study the sensitivity of our results to the number of constituents in each proton in Sec. V. Finally, our conclusions and future lines of work are given in Sec. VI.

II. WOUNDED HOT SPOTS MODEL
In this section, we briefly expose the main ingredients of our Monte Carlo Glauber event generator for protonproton interactions that follows similar steps than others in the literature [30,31]. For a detailed description of the model see [1]. In the following, a proton is considered to be formed by 3 hot spots. The comparison of our results for N hs = (2, 4) are given in Sec. V.
In each p+p event, after generating a random impact parameter for the collision, we sample the transverse positions of the three hot spots in each proton { s i } according to the distribution where C is a normalization constant and R is the aver-age radius. It should be noted that the extension of Eq. 3 to an arbitrary number of hot spots is direct. Most of the models in the literature [12,20,23,[32][33][34][35] implement a proton geometry following the two first terms of Eq. 3 i.e. the hot spots are distributed according to Gaussian functions with the natural constraint of fixing the centre of mass of the constituents system to the centre of the proton. However, with this set up the most probable configurations are the ones with three hot spots in the middle of the proton and the one with two hot spots fully or partially overlapping and the third one separated due to the δ-function like in the quark-diquark model. The third term of Eq. 3 allows us to go beyond these approaches by implementing short range repulsive correlations among all pairs of hot spots that effectively enlarge the mean transverse separation | s i − s j | between them. The size of this correlation is controlled by an effective repulsive core distance r 2 c ≡ R 2 /µ. The original motivation to consider these additional spatial correlations within our model was their critical impact on the dynamical explanation of the hollowness effect [25]. Although it constitutes the main novelty of our phenomenological model with respect to others in the literature in this context it should be noted that the necessity of spatial correlations has been already entertained in the nuclear case [36][37][38]. Further, a similar mechanism prevents the ropes in the DIPSY event generator to be in a highly energetic color state [39]. All along the manuscript we will focus on comparing the results obtained in the uncorrelated scenario (µ → ∞) with the correlated case.
Once the hot spots are located in both target and projectile, the next step in the Monte Carlo simulation is to decide which of them have been wounded [40,41] i.e. have collided at least once. Our collision criterion consists of sampling the inelasticity density that depends on the radius of the hot spot R hs , the transverse distance between the pair of hot spots considered d and the ratio of real and imaginary parts of the hot spot-hot spot scattering amplitude ρ hs . Thus, in each event, the maximum number of wounded hot spots N w and collisions N coll is 6 and 9 respectively. In Section IV A we will use these two variables to characterize interaction topologies. Subsequently, we consider that each wounded hot spot located at (x w , y w ) deposits a random amount of entropy according to where s 0 fluctuates independently for each wounded hot spot following a double Gamma distribution where parameters, given in Table II, are fixed by the assump-tion that entropy deposition is related to particle production [1]. Moreover, the centrality classes considered At this point, all the wounded hot spots contribute to the calculation of the spatial eccentricity moments that characterize the initial geometry anisotropy of the collision where · in Eq. 6 denotes the average weighted by the entropy deposition given by Eq. 5. We compute Eq. 6 in the participant plane on an event-by-event basis.
Regarding the model parameters, by default we use the same values as the ones from [1] at √ s = 13 TeV given in Table I. However, in Section IV B we extend our calculation to other regions of the parameter space. The three correlation scenarios under consideration in this work are the following: first, r c = 0.4 refers to the correlated scenario with {R hs , R, ρ hs } constrained to reproduce the extrapolated values of the total cross section and the ratio of real and imaginary parts of the scattering amplitude [42]. In the second case we impose the latter constraints to {R hs , R, ρ hs } after setting r c = 0. However, in order to do a more realistic comparison between the correlated and uncorrelated scenarios it is necessary to avoid the intrinsic swelling effects due to the presence of repulsive correlations. The case labeled as " s 1 fixed" constitutes an attempt to perform this task by fixing the r.m.s of the spatial probability distribution given by Eq. 3 to be the same as in the r c = 0.4 scenario. In the following plots the uncorrelated results will be exhibited as a band bounded by the r c = 0 and s 1 fixed cases to display the different possibilities considered.
Once the building blocks of our model have been presented in the next sections we display its results for the normalized symmetric cumulants given by Eq. 2 obtained after generating 4.5 million events. Only the events with at least two wounded hot spots contribute to the averages in the following plots.

III. NORMALIZED SYMMETRIC CUMULANTS VS CENTRALITY
The most important result of this paper is shown in Fig. 1 where we represent the event-averaged value of NSC(2,3) as a function of centrality. A common feature in the three correlation scenarios is the fact that   NSC(2,3) decreases from peripheral to central collisions as suggested by data. Focusing on the effect of the shortrange repulsive correlations we observe how they enlarge the positive correlation of ε 2 and ε 3 in the peripheral regime. However, their repercussion in the very central collisions is precisely the opposite. Finally, the most striking effect of the spatial correlations is observed in the ultra-central bins [0-0.1%] and [0.1-1%]: only in the r c = 0.4 case there exists an anti-correlation of ε 2 and ε 3 as data dictates. Then, we conclude that the experimental evidence of NSC(2,3)< 0 may back up the necessity to consider correlated proton constituents. An important comment is in order at this point. A direct comparison with the experimental data is not straightforward especially in the low multiplicity regime where non-flow dijet contributions, totally absent in our initial state coordinate space approach, dominate the measured values of NSC(n,m). Another important issue is to ensure that the centrality bin selection is exactly the same both in our approach based on the entropy deposition and in the data in terms of N offline trk [7]. Extending our calculation to higher centrality bins is doable but computationally expensive. However, a precise theory-to-data comparison, although desirable, is not the main objective of this work. Our purpose is to present, for the first time in the literature, a particular mechanism i.e. the presence of spatial correlations inside the proton that builds up a negative value of NSC(2,3) in the highest centrality bin at the geometric level.
In the case of NSC (2,4), the role of the repulsive correlations is qualitatively the same as in the NSC(2,3) calculation: in peripheral collisions the value of NSC (2,4) is larger in the r c = 0.4 case than in the uncorrelated scenarios and the situation gets reversed at barely the same centrality bin. As well, we find the absolute value of NSC (2,4) to be larger than NSC (2,3) in all the centrality bins as it is the case in the data. We would also like to remark that in our approach the symmetric cumulants are almost flat in the mid-to-peripheral interactions but thanks to a dissection of the very central bins we see a clear centrality dependence. This is consistent with our previous calculations of the average values of the spatial eccentricity moments [1].
A geometric and intuitive interpretation of the fact that only in the correlated case NSC(2,3)< 0 in the [0-1%] centrality bin is given in the following section. It should be noted that, for this purpose, we have merged the two highest centrality bins, [0-0.1%] and [0.1-1%], into a single one in order to improve the statistics.

IV. ULTRA-CENTRAL EVENTS
All the results presented in this Section refer to the [0-1%] centrality bin. We restrict our calculations to this bin because as we have emphasized in the previous section, we are interested in the change of sign of NSC(2,3).

A. Role of the interaction topology
In order to capture the effect of the spatial correlations we characterize each proton-proton interaction by its number of wounded hot spots and the number of collisions (N w , N coll ), the two basic quantities of any Monte-Carlo Glauber calculation. We dub each (N w , N coll )configuration as interaction topology. In our case, given that we consider the proton to be formed by three hot spots N w ∈ [2,6] and N coll ∈ [1,9].
We begin our analysis by computing the average number of collisions as a function of the number of wounded hot spots for the three different scenarios introduced above. The results are shown in Fig. 3. First of all, as we describe the entropy deposition in an incoherent way i.e. on average the more wounded hot spots the more entropy is deposited, the configurations in which only two hot spots collide cannot create enough entropy to be part of the [0-1%] centrality bin. Then, the minimum number of wounded hot spots is three and, in this case, N coll = 2 in all the correlation scenarios as it is the only existing configuration. However, for N w > 3 the average number of collisions starts to differ between the three different cases. We observe that N coll is systematically reduced when including repulsive correlations with respect to the uncorrelated cases. This effect has a very straightforward interpretation: enlarging the mean transverse distance between the hot spots reduces the probability of having interaction topologies with a high number of collisions. In other words, the repulsive correlations spread the hot spots in the transverse plane and as a consequence enhance the probability of the hot spots to collide by pairs over the configurations in which all hot spots in the projectile interact with all the others in the target, as it is schematically represented in Fig. 4.
To connect this fact with the total value of NSC(2,3) we would like to understand the individual contributions from the different interaction topologies. For this purpose we define a weighted version of NSC(n,m) denoted NSC w (n,m) as follows  NSC w (n, m) ≡ P(N w ) · P(N w |N coll ) · NSC(n, m) where • P(N w ) is the probability of having a certain number of wounded hot spots.
• For a given N w , P(N w |N coll ) represents the probability of having a certain number of collisions between the hot spots.
• NSC(n, m) Nw,N coll is the value of NSC(n,m) for each interaction topology.
The error of NSC w (n,m)is computed by adding the statistical uncertainties of each term in Eq 7 in quadrature. Essentially, by summing NSC w (n,m) over all the possible configurations (N w , N coll ) one recovers NSC(n,m). This new quantity allows us to decompose the value of NSC(2,3) and investigate the contribution of each interaction topology separately. From now on, to facilitate the discussion, we only show the comparison between s 1 fixed and r c = 0.4 scenarios. We have checked that the same conclusions as in the s 1 fixed case hold for r c = 0.
In Fig. 5 we show a particular example of the output of our calculation for NSC w (2,3) by selecting the events with N w = 6. Two important results can be extracted from this figure. First, as already suggested by Fig. 3, configurations with a large number of collisions, e.g. N coll > 6, only occur in the uncorrelated case where the three hot spots are closer to each other or, equivalently, clustered. Second, and more important, the value of NSC w (2,3) shows a clear dependence on N coll : configurations with a smaller number of collisions reduce the value of NSC(2,3) and, eventually, contribute negatively. Then, in our picture, the inclusion of spatial correlations inside the proton modifies the weight of each interaction topology in such a way that these configurations are enhanced. This feature provides a natural explanation for the different sign of NSC (2,3) in the uncorrelated and correlated scenarios. Far from being a casual coincidence or an artifact this effect is observed for any number of wounded hot spots as it is depicted in Fig. 6. In the top pannel we show the event-averaged value of NSC w (2,3) with respect to the number of collisions for N w = 3 to N w = 6 for the r c = 0.4 case. Once again, the configurations that contribute more to the total value of NSC (2,3) are the ones with a large number of wounded hot spots that interact a small amount of times. In opposition, as displayed in the bottom pannel, the interaction topologies that have associated a negative NSC(2,3) are extremely suppressed in the uncorrelated scenario where the configuration with the biggest weight and precisely positive value of NSC(2,3) is (N w = 4, N coll = 3).
Then, by computing NSC w (2,3) for the different interaction topologies we find that the origin of the negative sign of NSC (2,3) in the r c = 0.4 scenario is due to the

B. Scan of the parameter space
To conclude our study we check the sensitivity of the obtained results on the values of the model parameters. Thus, we focus on the correlated scenario and study the dependence of NSC(2,3) on the radius of the hot spot and the repulsive core distance in the [0-1%] centrality bin. As it could be argued that r c = 0.4 fm is a large repulsive distance that may be unrealistic we explore the results of our model for r c = 0.25 fm. In the case of R hs , we choose 4 different values in our scan {0.15, 0.25, 0.32, 0.4} fm. Consequently, the parameters of the Gamma distribution for the entropy deposition (see Eq.9 in [1]) are extracted in all the cases by fitting the experimental charged-particle multiplicity distributions P(N ch ) [43]. The other two parameters of our model, namely R and ρ hs , remain fixed to their default values given in Table I. Except in the chosen values for R hs and r c appearing in Table I, i.e. R hs = 0.32 fm and r c = 0.4fm, the requirement that our model reproduces the p+p total cross section and ρ is not fulfilled. Removing these phenomenological constraints allows to pinpoint the effect of just varying the radius of the hot spot or the correlation distance in our results. In Fig. 7 we represent the event-averaged value of NSC(2,3) as a function of R hs for the two different values of the correlation distance considered. First, we observe that by reducing the value of r c for a given value of R hs we get closer to the uncorrelated case and thus the value of NSC(2,3) is enlarged and pushed to the positive regime, as expected. However, this statement is not universal as it breaks down when R hs 0.22 fm. In this scenario of very small values of the radius of the hot spot, i.e. R hs 0.22 fm, ε 2 and ε 3 are positively correlated for both values of the repulsive core distance and the value of NSC(2,3) is larger in the r c = 0.4 case. This result indicates that NSC(2,3) is not sensitive to R hs and r c independently but to the interplay of both scales. In other words, NSC(2,3) depends on a generic function of the radius of the hot spot and the repulsive core distance f (R hs , r c ).
As a first and simple guess to the functional form of f (R hs , r c ) we choose it to be the ratio of the two scales involved i.e. f (R hs , r c ) = R hs /r c . This ratio has a transparent interpretation by characterizing the degree of repulsion: if R hs /r c 1 the shape of the proton resembles the uncorrelated scenario where the hot spots can largely overlap in transverse space. The results of NSC (2,3) for different values of R hs /r c are displayed in Fig. 8. We can distinguish three regimes in this plot. On the one hand, when R hs /r c ≥ 1 the geometric picture of the proton approaches the uncorrelated scenario and the value of NSC(2,3) increases monotonically starting to be positive for R hs /r c 1.3. Moreover, when R hs /r c = 1 the value of NSC(2,3) is identical within error bars for both correlation scenarios r c = 0.25 and r c = 0.4 supporting the idea that NSC(2,3) depends on f (R hs , r c ) = R hs /r c . However, on the second regime characterized by 0.6 R hs /r c 1, we find an abrupt change of the value of NSC(2,3) when slightly increasing the ratio R hs /r c from 0.6 to 0.63. This suggests that a residual dependence of NSC(2,3) on the other scale of the problem, R, may exist. Finally, configurations in which the hot spots are much smaller than the repulsive core distance between them i.e. R hs /r c 0.6 result into a positive correlation between ε 2 and ε 3 . Then, our study favors values of 0.6 R hs /r c 1.3 in order to be compatible with the experimental observation of NSC(2,3)< 0 in the highest centrality bin. Unfortunately, this interval is large enough to be compatible with a picture of the proton in which the hot spots transverse separation is larger than in the uncorrelated case but still they can overlap (R hs /r c ∼ 1.3) and with a much more dilute description in which the probability of two hot spots to overlap is highly supressed (R hs /r c ∼ 0.6).

V. SENSITIVITY OF NSC(2,3) TO N hs
All along the manuscript we have considered that the proton is constituted by 3 gluonic hot spots. This is the most natural scenario when a direct correspondance between the Fock space of valence partons and the hot spots is assumed. However, this relation is arguable as, while being extensively used as a phenomenological tool, the ultimate dynamical origin of the hot spots remains as an open debate. Therefore, it is opportune to check the reliability of our results after variations of this   parameter, N hs .
Following the ideas of the previous sections, we focus our discussion on the results for NSC(2,3) after considering the two more straight-forward extensions of our model: N hs = 2 and N hs = 4. In order to make a fair comparison between the three different scenarios, N hs = (2, 3, 4), we choose representative values of the parameters {R hs , R} that fulfill two constraints. As in the previous sections, the experimental value of the total p+p cross section is reproduced. Further, the proton radius defined as R p = √ N hs s 2 1 + R 2 hs , where s 1 is the r.m.s of the spatial probability distribution given by Eq. 3, should not depend on the number of hot spots that the proton contains so we fix it to be the same in all the cases. All in all, the values of the parameters for the correlated and s 1 fixed cases are given in Tables III and  IV. We start by exploring the dependence of the eventaveraged value of NSC(2,3) on the number of hot spots in the correlated scenarios as displayed in Fig. 9 as a function of centrality. The differences between the three cases start to appear in mid-to-ultra central collisions. There exists a clear trend towards smaller values of NSC(2,3) when a bigger number of hot spots is considered. Specifically, the negative sign of NSC (2,3) in the high centrality bins is not achieved when N hs = 2 even with correlations. Thus, we conclude that with the selected parameters the minimum number of hot spots to describe the onset of the anti-correlation between ε 2 and ε 3 is N hs = 3. It should also be noted that the inclusion of an additional hot spot i.e. N hs = 4 helps to make NSC(2,3) even more negative in the highest centrality bins although the effect is small when compared to the drastic impact of changing from N hs = 2 to N hs = 3.
To conclude the discussion on the sensitivity of our model to the number of gluonic hot spots, in Fig. 10 the comparison between the correlated and s 1 fixed scenarios for N hs = 2 (top) and N hs = 4 (bottom) is displayed. First and foremost, the effect of including spatial correlations is invariant under changes in the number of hot spots: in peripheral collisions they enlarge the positive correlation of ε 2 and ε 3 while favoring a negative sign of NSC(2,3) with respect to the uncorrelated scenario. Albeit the correlated curve is always below the uncorrelated scenario in the highest centrality bins an important comment is in order: NSC(2,3) is compatible with negative values, within statistical uncertainty, in the [0-0.1%] bin for the uncorrelated case. This fact reinforce the idea remarked in Sec. IV B: the interplay of the different scales {R hs , r c , N hs } is decisive in the sign of NSC(2,3) within our framework. For N hs = 4 the weight of the configurations with a large number of wounded hot spots and a small number of collisions is large enough so that the spatial correlations are not essential to obtain a negative NSC (2,3) in the [0−0.1%] bin. However, as these configurations are enhanced in the correlated scenario, the anti-correlation of ε 2 and ε 3 is stronger than in the uncorrelated case just as in the case of N hs = 3. Finally, we have checked that the effect of spatial correlations on the average values of ε 2 and ε 3 is qualitatively the same as in the N hs = 3 case studied in [1]. To sum up, although the negative sign of NSC (2,3) in the highest centrality bins is not a unique feature of the correlated scenario but relies on the interplay of the different scales, the inclusion of repulsive correlations provides a mechanism to reduce its value in the highest centrality bins.

VI. DISCUSSION AND OUTLOOK
The experimental finding of a resembling behavior of the fluctuations of the Fourier harmonic coefficients v n , in terms of the symmetric cumulants SC(n,m) in p+p, p+Pb and Pb+Pb collisions at LHC energies constitutes a new piece of the puzzle on whether collective effects are being observed in small systems. Data suggests that the fluctuations of v 2 and v 3 are anti-correlated i.e. SC(2,3)< 0 in proton-proton collisions with N offline trk ∼ 100. In this very-high multiplicity regime, where the non-flow contributions arising from jet correlations are subleading, the value of the symmetric cumulant is sensitive to the initial state fluctuations of the collision.
In this article, we perform a systematic study of the normalized symmetric cumulants NSC(2,3) and NSC (2,4) in proton-proton interactions at √ s = 13 TeV in terms of the initial geometry of the collision. We rely on the wounded hot spot model used in [1] where the proton is regarded as a system of three gluonic hot spots whose transverse positions are not independent but correlated via a repulsive core distance. We find that the effect of these repulsive correlations is key, within our relatively simple geometric model, to obtain an anti-correlation of the fluctuations of ε 2 and ε 3 in the highest centrality bin. More precisely, they enhance the probability of having interactions with a large number of wounded hot spots colliding a small amount of times. These interaction topologies are responsible in our set up of the negative sign of the symmetric cumulant NSC (2,3) in the correlated scenario. Further we explore the dependence of NSC(2,3) on the values of the repulsive distance and the radius of the hot spot concluding that it is not a function of these two variables independently but to a combination of both such as their ratio R hs /r c . We also show that the values of NSC(2,3) are sensitive to the variation of the number of hot spots that constitute the proton. Specifically, we find that, within the correlated scenario, adding an extra hot spot to our description reinforces the negative sign of NSC(2,3) while reducing it to N hs = 2 pushes NSC(2,3) towards positive values. Moreover, the enhanced probability of having configurations with a large N w /N coll when N hs = 4 permits a negative value of NSC(2,3) within the uncorrelated scenario although its absolute value is smaller than in the correlated case. This fact reinforces the argument that the sign of NSC(2,3) is sensitive to the interplay of the different scales.
Our study confirms the idea that NSC(2,3) in protonproton interactions is extremely sensitive to the initial state fluctuations. Further, as we have shown, it can help to discriminate between different parameterizations of the proton's geometry and, concretely, to realize the importance of spatial correlations. A precise characterization of the proton's geometry has a direct impact in the flow studies in proton-proton collisions. Furthermore, parametrizing how are the subnucleonic degrees of freedom arranged inside the proton is an essential input in event generators that aim to describe multi-parton interactions, the mechanism that dominates the underlying event at LHC energies.
In order to confirm the conclusions exposed in this manuscript the natural continuation of this work would consist on feeding a relativistic viscous hydrodynamic simulation with our initial entropy density profiles and check if the effect of these spatial correlations is washed out by the evolution or, on the contrary, it impacts the values of the Fourier flow coefficients.