Signals for fluctuating constituent numbers in small systems

TL;DR

This paper extends the TRENTo model to include fluctuating numbers of partons in nucleons, incorporating PDFs and momentum transfer dependence, to better understand initial conditions in heavy-ion collisions.

Contribution

It introduces a novel sampling method for parton numbers based on PDFs within the TRENTo framework, accounting for momentum transfer effects.

Findings

Distributions depend on momentum transfer and scale.

Observable sensitivities to parton number fluctuations are analyzed.

The model captures fluctuations in initial conditions more realistically.

Abstract

We propose an extension of the initial condition model TRENTo for sampling the number of partons inside the nucleons that participate in a heavy-ion collision. This sampling method is based on parton distribution functions (PDFs) and therefore has a natural dependence on the momentum transferred in the collision and the scale being probed during the collision. We examine the resulting distributions and their dependence on the momentum transfer. Additionally, we explore the sensitivity of different observables on the number of partons using the TRENTo framework and the estimators available therein for final-state observables.

Equations

$$Q^2=sxy$$

$$x\ge \frac{Q^2}{s}=x_{min}.$$

$$N_{\text{ch}}$$

$$⟨p_{\mathrm{T}}⟩$$

$$v_n$$

Full text

I Introduction

One of the major tools for enhancing our understanding of quantum chromo dynamics (QCD) are ultra-relativistic heavy-ion collisions. These collisions, conducted mainly at the large hadron collider (LHC) at CERN and the relativistic heavy-ion collider (RHIC) at Brookhaven National Laboratory Harris:1996zx; PHENIX:2003iij; PHENIX:2003qra; STAR:2008med; ALICE:2013mez; ALICE:2022wpn; Harris:2023tti, allow to probe the properties of hot and dense QCD matter by comparison between experimental data and theoretical models Bernhard:2016tnd; Bernhard:2019bmu; Nijs:2020roc; JETSCAPE:2020mzn. A majority of these models use a fluid-dynamic description for the bulk part of the collision and encode the underlying QCD dynamics through the equation of state and transport coefficients, such as shear and bulk viscosities. While it is possible to obtain the equation of state and the transport coefficients through means such as lattice QCD, renormalization group equations and others HotQCD:2014kol; Christiansen:2014ypa; Policastro:2001yc, the description of the first few moments of the collision from first principles is much harder. Despite recent advances in this area Schenke:2012wb; Niemi:2015qia; Kurkela:2018vqr, the precise dynamics is not yet fully understood and computationally expensive to study. A popular alternative choice is the usage of initial state models based on more parametric ansatzes, such as Monte Carlo Glauber models Miller:2007ri; Loizides:2017ack or the T${}_{\text{R}}$ENTo initial condition model Moreland:2014oya; Ke:2016jrd; Moreland:2018gsh; Soeder:2023vdn. Since these models include parameters that may be difficult to directly relate to calculable QCD quantities, a common way of determining their value is by calibrating the respective models to experimental data. Depending on the type of measurement at hand, a different set of parameters can be constrained based on the model-data comparison.

One recent example regarding the discovery of new parameters for calibration was the search for the chiral magnetic event at RHIC STAR:2021mii through the collision of two isobar isotopes (${}^{96}_{44}\mathrm{Ru}$ and ${}^{96}_{40}\mathrm{Zr}$). Extensive studies have demonstrated a strong dependence of the ratios between isobars on nuclear structure parameters, such as the neutron skin thickness and deformations Giacalone:2025vxa; Xu:2021uar; Xu:2021vpn; Jia:2021oyt; Nijs:2021kvn; Li:2018oec.

While the the modeling of deformed nuclei is well established and the relation between the nuclear deformation and their effects on flow coefficients in particular has been demonstrated Giacalone:2024luz, the modeling of sub-nucleonic degrees of freedom for hydrodynamic initial conditions remains an active area of research in the heavy-ion community Moreland:2018gsh; Nijs:2020roc; Soeder:2023vdn; Shen:2017bsr; Giacalone:2021clp. A common way for including sub-nucleonic structure in the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo model is the addition of $m$ partons inside each proton and neutron with $m$ ranging from three to values of around ten, with $m$ being a fixed number for all nucleons in both ions. However, from deep inelastic scattering and other experiments, it is known that protons and neutrons have a rich substructure, especially when probing them at small Bjorken $x$, also known as saturation region Kovchegov:2014kua. This dynamics is encoded in so-called parton distribution functions Lai:1999wy, which give the probability to find a parton at a given Bjorken $x$. Additionally, the number of partons for any nucleus at a given momentum transfer is not fixed, but fluctuates on an event-by-event basis. In this paper we present a novel way of sampling the number of constituents of a nucleon in T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo , based on parton distribution functions, accounting for event-by-event fluctuations of the number of partons.

The sampling algorithm together with the resulting distributions are discussed in Sec. II. The signature of these fluctuations is then examined and shown to be strongest in small collision systems, such as oxygen or neon, (Sec. III). Additionally, we examine different observables and determine a combination of known observables displaying an enhanced sensitivity to the parton number in Sec. III.

II Parton number sampling algorithm

Before discussing the details of the sampling algorithm, we will first discuss how sub-nucleonic structure is currently treated in T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo : In T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo the positions of the nucleons inside the nuclei are sampled based on the nuclear density. The probability for a collision is then determined based on the cross section and distance for each nucleon pair in the two colliding nuclei. Since this collision process is so far not well understood from first principles, the amount of energy deposited in each nucleon-nucleon collision is modeled by placing a Gaussian at the collision site, whose normalization is being drawn from a Gamma distribution to allow for fluctuations. At this step the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo model allows for the inclusion of subnucleonic degrees of freedom, i.e. quarks and gluons, by means of depositing not one large Gaussian at the collision site, but rather depositing $m$ smaller Gaussians around it 111The sampling for their positions is designed to yield the same nucleon size as the single Gaussian in the limit of large $m$. It was already demonstrated that this fixed number of sub-nucleonic degrees of freedom can be constrained using bayesian inference with $m$ found to be between four and five222The number of constituents is to be understood as rounded value of the reported maximum likelyhood value of the posterior. with rather large uncertainties Nijs:2021clz; Nijs:2020ors; Nijs:2020roc. Additional, stronger constraints can arise in smaller collision systems Huang:2025cjm.

This model of a fixed number of hotspots at each collision site only superficially reflects our understanding of the structure of the proton and neutron for two main reasons:

Firstly, it is well established that protons and neutrons have complicated and rich substructure consisting out of three valence quarks together with a fluctuating amount of sea quarks and gluons, resulting in a fluctuating parton number when probing the same proton or neutron multiple times. The second discrepancy between the model and our understanding is the fact that if the proton had a fixed number of partons, not all would need to interact in the inelastic collision due to the different energy (and therefore time) scales set by the inelastic scattering compared to the scales in parton-parton interactions. Both of these effects lead to a fluctuating number of partons in each nucleon-nucleon collision.

To determine this fluctuating number of partons, we propose a novel distribution based on the parton distribution functions (PDFs) describing the momentum distribution of partons inside of the nucleus. The parton distribution functions hereby serve as probability distributions for finding a parton of a certain type in the momentum range from $x$ until $x+\mathrm{d}x$, with $x$ being the momentum fraction carried by the parton of the total momentum of the proton or neutron. Note that this is not a probability distribution function in the strict sense, since it is not normalized to unitiy for the different species. Nevertheless, we will use it to obtain the parton number for a given event based on the moment fractions sampled from the PDFs:

For the sampling of the parton number, we will assume that the three valence quarks found in the proton and neutron are strongly correlated among each other and therefore always participate in the reaction. Therefore, the sampling will begin by drawing the momenta of the three valence quarks from their respective up and down quark distribution functions. The sampled momenta are then added up, $\sum_{val}x_{val}=x_{tot}$. Afterward, depending on the momentum transfer $Q^{2}$ the species for the next sample is determined randomly from gluons, up, down, strange, charm and bottom quarks, where the heavier species are only sampled for large enough $Q^{2}$. Subsequently, either the momentum of one gluon or a quark anti-quark pair is sampled from the respective distribution. The drawn momenta are then added to the total momentum sum $x_{tot}$. This procedure beginning from the random species determination is then repeated until the total momentum sum exceeds one: $x_{tot}>1$. The total number of partons is then given by three for the valence quarks plus the amount of drawn sea gluons, quarks and anti-quarks. At this point it is important to note that the number of partons will be bounded from below, despite PDFs sharply rising for $x\to 0$, due to a lower cut-off imposed on $x$ by kinematic considerations:

Using the relation between the momentum transfer $Q^{2}$, the center of mass energy $s$, the target parton momentum fraction $x$ and the projectile energy loss $y$

[TABLE]

and using $y\leq 1$, we find the minimal value for $x$, given by

[TABLE]

In spite of the small value of $x_{min}\approx 10^{-6}$ for top LHC energies, it nevertheless provides a cut-off for the parton distribution function, ensuring a finite maximal value, see Fig. 1. While the parton distribution functions sharply rise for $x\to 0$, the cutoff $x_{min}$ ensures a finite maximal values for each distribution, resulting in a normalizable distribution. Furthermore, since the cut-off point is larger than zero $x_{min}>0$, it ensures a maximal value for the number of partons given by $m_{max}=1/x_{min}$. The parton distribution functions used for the sampling are taken from the UCL HEP group Harland-Lang:2014zoa 333We want to highlight that ideally, the partons would be sampled using the $n$-parton distribution function (instead of the 1-parton distributions), including all correlations between the partons. However, since these are not yet known, we will sample partons individually, neglecting some of their correlations..

The resulting distribution for the parton number $m$ for different values of $Q^{2}$ are displayed in Fig. 2. All distributions are bounded from below by $m=3$ through the presence of high momenta valence quarks. Since the PDFs are less peaked for lower values of $Q^{2}$, there are more high momenta valence quarks at low $Q^{2}$, resulting in a distinct peak of the parton distribution at $m=3$. Since a collision with a larger momentum transfer probes the system at smaller length scales, the number of parton grows and therefore the maximum of the parton distribution also shifts to larger values.

Similarly to the maximum of the distribution, also the average parton number grows with $Q^{2}$ initially, see Fig. 3. However, for large enough $Q^{2}$, the average number of partons starts to decrease again. This can be attributed to the interplay between the scale being resolved in the scattering process getting smaller with larger $Q^{2}$ and the momentum cutoff getting larger with $Q^{2}$. The resolution of smaller length scales leads to the PDFs being more peaked toward $x\to 0$, resulting in more low momentum partons being sampled, increasing the total number of partons. However, with rising $Q^{2}$, the excluded momentum range below $x_{min}$ also grows, leading effectively to partons with higher momenta, therefore reducing the parton number (see Fig. 5).

This interplay between the changing distribution functions and the cutoff with changing momentum transfer also leads to a decreasing standard deviation of the parton distributions with increasing $Q^{2}$, as shown in Fig. 4. The higher cutoff reduces the range the sampled momentum is able to fluctuate in, leading to a smaller standard deviation of the distribution of partons.

Since this method of obtaining the parton number solely relies on parton distribution functions, it is possible to employ this method for different collision systems and energies through the use of DGLAP equations Gribov:1972ri; Dokshitzer:1977sg; Altarelli:1977zs. An additional extension using the different distributions for protons and neutrons is straight forward by using the respective proton or neutron distribution functions. In a similar fashion, the difference between different nuclei can be accounted for by using the nuclear distribution function instead of the proton or neutron ones.

Since the parton number can get large, especially for large center of mass energies, a large number of sampling processes can be required, resulting in a higher evaluation cost of initial state model, especially for heavy nuclei. A possibility to circumvent this added cost to the evaluation of the initial conditions is to employ a parameterization or fit of the parton distributions presented in Fig. 2. A crude parametrization of the distribution can be obtained from its cumulants, with the first two (mean and standard deviation) being given in Fig. 3 and Fig. 4 as function of $Q^{2}$.

Additionally, the question arises, if one should take all of the sampled partons into account. Despite the large correlations between the individual partons in the nucleus, not all of them might participate in the interaction due to possibly very large virtualities. In the entirety of this work, we assume the correlations to be strong enough for all particles, no matter their virtuality, to interact. However, a further examination of such a ”cutoff” in virtuality presents an intriguing possibility for the extension of this work.

Another extension of this sampling method is the usage of so-called transverse momentum distribution functions (TMDs) instead of parton distribution functions. Transverse momentum distribution functions can be understood as extensions of parton distribution functions by including additional information about the transverse momentum distribution of the individual partons. This additional information can be used to sample the transverse momenta of the partons participating in the collision. The inclusion of this information can then allow for the calculation of a initial non-zero fluid velocity.

While we believe that this model has a strong physics motivation and propose its implementation in the next version of an initial state model, such as T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 444The actual implementation requires additional considerations regarding the amount of deposited energy per parton depending on their respective momentum fraction., we will use calculations at fixed $m$ to examine the sensitivity of different observables on the number of partons using initial state estimators.

III Sensitivity of different observables

In this section we want to gauge the sensitivity of different observables and collision systems on the proposed fluctuating number of constituents. In principle, every combination of observables can be used to constrain the number or distribution of partons inside each nucleon. However, some of these observables will show greater sensitivity than others. To examine these effects, we will make use of the well known proportionalities Giacalone:2020ymy linking initial state quantities, such as the entropy, energy and eccentricities to the number of charged particles, the mean transverse momentum and the harmonic flow coefficients,

[TABLE]

Even though these are not precise relations, they nevertheless can aid in gauging the sensitivity to different observables by employing only initial state simulations, which are relatively cheap compared to a full fluid dynamical simulation.

We already can establish tendencies of the sensitivity of different observables: Since the amount of energy carried by each nucleus is the same, irrespective of the number of partons, the energy deposited in each collision will remain similar. Likewise, the entropy of the collision will remain largely unchanged, although one could make the argument that the strong correlations between the individual partons will increase the entropy of the initial profile. Since the description of this entropy from first principles is non-trivial, we will assume that it can be neglected in this context. Having multiple constituents inside each nucleon might lead to more defined geometric features of the initial profile compared to the case of larger, smooth nuclei, resulting in a change of the initial eccentricities. Therefore, the largest signal is expected in the eccentricities.

Depending on the centrality and the collision system, up to 2000 binary nucleon-nucleon collisions take place in a single heavy ion collision, effectively sampling the distribution of constituents for each of the binary collisions. Regrettably, it is not possible to access these 2000 samples individually, which would allow to put precise constraints on the distribution and its moments. Since all of the individual energy and entropy deposition sites are combined into one initial energy/entropy density profile, even the initial condition used for fluid-dynamic simulations only contains an average of the partonic distribution function. More of this information will be washed out during the hydrodynamic evolution phase and subsequent particle production. Since this effective averaging scales with the number of binary collisions, we expect more insights into the parton distribution when considering events with a lower amount of binary collisions. This can either be achieved by considering more peripheral events (see Fig. 7) or by examining smaller collision systems. Interestingly, also ultra-central events show some dependence of the number of partons, even if only a slight one, as depicted in Fig. 6. This stems from the number of spectators going to zero, resulting in a similar value of $N_{\text{coll}}$ for each event equaling no $N_{\text{coll}}$ fluctuations. When using a fluctuating distribution for the number of partons, this would result in probing the average of the distribution.

The expected difference in the eccentricities is demonstrated in Fig. 8, where the eccentricities as function of centrality are displayed for three different systems: PbPb, OO and NeNe collisions. Similarly to the initial entropy small deviations for the very central collisions also appear in the eccentricities. These differences are then much more pronounced for peripheral collisions. Here it is interesting to note that the even harmonics, i.e. $\varepsilon_{2}$ and $\varepsilon_{4}$, show less separation, especially for the large centralities, compared to the odd harmonics $\varepsilon_{3}$ and $\varepsilon_{5}$. This different behavior at large centrality is a result of the even harmonics mainly being driven by the geometric shape of the overlap area, whereas the odd harmonics are mainly driven by fluctuations, which get enhanced through the added nuclear substructure. Since a higher number of binary collisions tends to wash out the details of the sub-nuclear structure, the difference for different $m$ appears at different centralities when considering a large system, such as a PbPb collision, compared to smaller systems, such as OO or NeNe collisions. These differences arise at around the same value of $N_{coll}$ in all systems, therefore appearing at smaller centralities for the smaller systems. To be able to study the effects of different $m$ we therefore will focus on OO collisions when examining the harmonic eccentricities further.

In order to enhance effects appearing due to differences in $m$ and suppressing common effects, we consider the ratios between the eccentricities for OO collisions, as depicted in Fig. 9555The differences between the various harmonics display a similar behavior.. All ratios, including even eccentricities, are again mainly driven by the geometric shape of the collision and therefore only show differences for different $m$ at large centralities. However, the ratio of the odd eccentricities is driven by fluctuations and therefore displays differences already around $50\%$ centrality.

The ratio in odd harmonics shows additional features, which could help in constraining the number of partons: While being a relatively flat curve for small values of $m$, the curve becomes steeper for larger values of $m$, with its maximum and minimum moving apart. Therefore, the smaller ratios at lower centralities could be leveraged together with the slope when going to larger centralities to constrain the value of $m$ or its distribution. This behavior can be attributed to $\varepsilon_{3}$ and $\varepsilon_{5}$ being relatively small for the central collisions, with the triangularity becoming more pronounced at higher centralities due to the sharper profiles at large values of $m$.

This feature is present in both OO collisions and NeNe collisions (see Fig. 10) and therefore could be an interesting addition to the analysis of the current OO and NeNe runs at LHC and possible future light ion runs.

While similar features also appear in other combinations of the eccentricities, such as their ratios, it is sufficient to include one additional observable containing the first two odd harmonics. To reduce uncertainties on the experimental side, we propose the above discussed ration between $\varepsilon_{3}$ and $\varepsilon_{5}$. Here, it should be noted that the fifth flow harmonic is not directly proportional to $\varepsilon_{5}$ anymore, due to non-linear flow effects. Therefore, the observed signal in the flow coefficients might not be as strong in the eccentricities.

IV Conclusion

In this work, we proposed a novel way of obtaining the number of partons in each proton and neutron participating in a heavy-ion collision. This method is based on the information obtained, among other through DIS experiments in parton distribution functions and therefore encodes the underlying QCD dynamics governing the structure of the proton. We examined the resulting distribution as a function of momentum transfer and discussed possible implications for the collision, while ensuring a finite number of partons for each sampled proton and neutron. Additionally, we propose possible modifications of this sampling algorithm accounting for different virtualities of the sea-quark and their antiparticles. As final part of the discussion of the sampling method, we sketch the possibility of obtaining a non-zero fluid velocity profile using additional information encoded in transfer momentum distribution functions.

In the second part of this work, we examined the dependence of different observables on the number of partons $m$ while keeping $m$ fixed. This is achieved by using the initial state estimators for particle number, average transverse momentum, and flow coefficients, namely entropy, energy, and eccentricities. We demonstrate that most observables show little to no sensitivity to the number of constituents, as expected through the inherent averaging dynamics of a heavy-ion collision. However, we were able to demonstrate that the odd eccentricities show some sensitivity to the value of $m$ at small centralities and a comparatively large sensitivity at large centralities. This sensitivity can be increased by considering their ratio, such that a clear change for different values of $m$ can be established already at around $50\%$ centrality.

We propose the inclusion of the ratio of the odd eccentricities as an additional variable in the upcoming analysis of small collision systems, as well as the addition of a fluctuating parton number based on PDFs in the next generations of initial state models, such as T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo .

Data availability statement

The data used for this research is available at github.

Acknowledgement

The authors want to thank Jyotirmoy Roy and Marston Copeland for useful discussions. This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, grant No. DE-FG02-05ER41367.