Glauber Monte Carlo predictions for ultra-relativistic collisions with ${}^{16}{\rm O}$

We explore Glauber Monte Carlo predictions for the planned ultra-relativistic ${}^{16}{\rm O}$+${}^{16}{\rm O}$ and p+${}^{16}{\rm O}$ collisions, as well as for collisions of ${}^{16}{\rm O}$ on heavy targets. In particular, we present specific collective flow measures which are approximately independent on the hydrodynamic response of the system, such as the ratios of eccentricities obtained from cumulants with different numbers of particles, or correlations of ellipticity and triangularity described by the normalized symmetric cumulants. We use the state-of-the-art correlated nuclear distributions for ${}^{16}{\rm O}$ and compare the results to the uncorrelated case, finding moderate effects for the most central collisions. We also consider the wounded quark model, which turns out to yield similar results to the wounded nucleon model for the considered measures. The purpose of our study is to prepare some ground for the upcoming experimental proposals, as well as to provide input for possible more detailed dynamical studies with hydrodynamics or transport codes.


I. INTRODUCTION
In the continued quest [1] for deeper understanding of the rich physics unveiled by ultra-relativistic nuclear collisions, proposals have been made to study collisions with 16 O beams, both at the LHC (see Sec. 9.10 of [1] for 16 O+ 16 O and Sec. 11.3 or p+ 16 O, describing the experimental programs that could be carried out in runs beyond the year 2022) and at RHIC [2].
Investigations of 16 O+ 16 O are motivated by the exploration of emergence of collectivity in small systems, which has attracted a lot of attention over the past few years [3][4][5][6][7][8][9][10]. This is a major issue, as it concerns the very nature of the initial dynamics in the created fireball (for recent overviews see, e.g., [11][12][13] and references therein). A distinct feature of the 16 O+ 16 O collisions is that with a similar number of participants as in the earlier studied p+Pb collisions, the participants are distributed more sparsely in the transverse plane. This is expected to lead to different subsequent evolution.
Studies of p+ 16 O collisions find a broader justification from the physics of air showers generated with cosmic rays [14] and our lack of full understanding of the production process, e.g., the cosmic ray neutrino puzzle (see [15] and references therein). They also carry significance for investigating the onset of collectivity in ultra-relativistic nuclear collision.
The purpose of this work is to provide some model predictions for the planned reactions that could be used in preparatory analyses for the experimental proposals. We use the Glauber [16,17] modeling, which has become a basic tool to describe the initial state due to its simplicity and phenomenological success. Our simulations are carried out with GLISSANDO 3 [18].
We note that an analysis similar to ours has recently been carried out for the LHC energies by Sievert and * Maciej.Rybczynski@ujk.edu.pl † Wojciech.Broniowski@ifj.edu.pl Noronha-Hostler [19], where the T R ENTo code [20] has been used for the initial conditions and hydrodynamics run with v-USPhydro [21]. Studies based on AMPT model [22] were presented by Huang, Chen, Jia, and Li in [23] for the RHIC collision energies. The details of our model implementation concerning the distribution of nucleons in 16 O as well as the NN reaction features are different from the above-mentioned approaches. Consequently, also the studied eccentricity measures differ to some extent, providing an independent estimate for model uncertainties in physical predictions for the considered reaction.
As the further evolution of the system with hydrodynamics or transport is outside of our present scope, we focus on flow observables which are not strongly sensitive to the hydrodynamic response, such as ratios of various flow coefficients or the normalized symmetric cumulants.

II. STRUCTURE OF 16 O
Before we embark on collisions of 16 O, it is worth to focus on its nuclear structure. This is relevant, as properties of collisions reflect the features of the projectiles (as well as, of course, the NN collision mechanism). Needless to say, the size of the nuclei affect the total cross section, whereas two-body correlations influence to some degree the flow observables [24]. To have the possibly most realistic 16 O nucleus, rather than using parameterizations of its one-body distribution, we take configurations from state-of-the-art dynamical nuclear physics calculations. Specifically, we use 6000 configurations from cluster Variational Monte Carlo (CVMC) simulations [25] with the Argonne v18 two-nucleon and Urbana X three-nucleon potentials, as provided in files in [26]. We stress that these dynamically generated distributions contain realistic nuclear correlations, which are absent (or put in by hand in the form of a repulsive core for the two-body distributions) when some simple parameterizations of the nuclear one-body distributions are used. The one-body nuclear radial densities are given in Fig. 1, where we plot the distribution of the centers of nucleons, ρ(r), conventionally in units of the central density ρ(0). The corresponding ms radius is r 2 = (2.6 fm) 2 , which folded with the proton charge form factor with ms radius r 2 p = (0.84 fm) 2 yields the ms charge radius of 16 O of which is comfortably in the right experimental range of (2.699(5) fm) 2 [27]. A standard measure of the nuclear two-body correlations is provided by the ratio between of normalized twobody probability distribution of nucleons in their relative distance r and the folding of their one-body distributions, . (2) In Monte Carlo simulations, this ratio is easily obtained by generating histograms in the relative distance between all nucleon pairs from the same configuration, divided by the "mixed" histogram, where the nucleons in the pair come from different configurations. By construction, the "mixed" two-body distribution exhibits no correlations and plays the role of the denominator in Eq. (2). The result of this procedure, obtained with GLISSANDO 3 [18], is shown in Fig. 2. We note a soft-core behavior at low separations r, where C(r) > 0 indicates repulsion. At larger r the correlations disappear, as expected. The noise in the figure is caused by rather low available statistics (6000 configurations from [26]).
In some cases, we will show comparisons to the case where the correlations are removed by the mixing method. The procedure used here is as follows: we take a nucleus whose nucleon positions are represented with spherical coordinates and regenerate randomly the angular coordinates, while retaining the radius. That way the correlations are removed, whereas the radial density distributions are preserved. This is important, as we do not wish to change size of the projectiles, which directly relates to the collision cross section. An equivalent method, yielding essentially the same results, is to form "mixed" nuclei with each nucleon taken from a different "physical" nucleus.
To summarize this section, in the following analysis of nuclear collisions we are going to use realistic dynamically generated configurations of 16 O, reproducing the charge radius and involving proper two-body correlations.

III. GLAUBER MODELING
The applied Glauber Monte Carlo approach (for a review see [28]) is described in detail in GLISSANDO 3 [18], so we are very brief here. Importantly, the adopted NN inelasticity profile, a.k.a. the Van Hove function [29,30], is obtained from fits of the COMPETE model parametrization implemented in the Particle Data Group review [31]. This parametrization provides the best available description of the pp and pp scattering over the range of the collision energies from √ s N N 5 GeV up to the highest LHC energies.
We use the wounded nucleon [32,33] variant of the Glauber model with an admixture of the binary collisions [34], which has been found necessary to describe the multiplicity distributions of the produced hadrons [35]. Thus the initial entropy deposition in the transverse plane, S(x, y), is proportional to a combination of the wounded and binary contributions controlled by the parameter α, where the densities ρ W (x, y) and ρ bin (x, y) are obtained from the positions of the point-like sources generated with the Glauber Monte Carlo, which are then smeared with Gaussian profiles of width 0.4 fm [18]. As is well known, the smearing effect quenches the eccentricities of the fireball. For the mixing parameter we take typical values, with α = 0.12 at √ s N N = 10 GeV, and α = 0.15 at √ s N N = 7 TeV and 10 TeV.
We note that the statistics of 6000 16 O configurations allows us to construct ∼18 M different collision events (not counting the random rotation of the nuclei and the change of the impact parameter), which is statistically more than sufficient for our studies. In the following we consider two collision energies for 16 O+ 16 O: √ s N N = 10 GeV, which is accessible in the beam energy scan at RHIC or at SPS, and √ s N N = 7 TeV, which has been studied at the LHC. The difference between these energies is in the value of the NN inelastic cross section (which grows from 31 mb up to 71 mb between these energies) and in the NN inelasticity profile [18].
First, we discuss the distribution of the number of wounded nucleons, N w , in "minimum bias" events, that is, for random unconstrained values of the 16 O+ 16 O impact parameter. The results for the two collision energies are compared in Fig. 3. We use the logarithmic scale, as typically done in experimental analyses. We note that at the higher energy the distribution is more flat at higher N w (with the obvious limit at N w = 32), since the higher value of the inelastic cross section makes it easier to wound more nucleons. The vertical lines indicate the corresponding centralities, obtained as quantiles of the distribution of N w .
This section contains the key results in view of the considered future 16 O + 16 O experiments. As this work is based of investigation of the initial condition obtained from the Glauber approach, we are going to focus on observables which are to a large degree independent of the hydrodynamic or transport expansion. This methodology is based on the shape-flow transmutation feature, appearing in hydrodynamic [36][37][38] or transport simulations [39], whereby the deformation of the initial transverse shape of the fireball leads to harmonic flow of the hadrons emitted at the end of the evolution. Moreover, the effect is manifest in an approximate proportionality of the flow coefficients v n to the initial eccentricities n , holding for n = 2 and 3 and for sufficiently central collisions: v n κ n n , (n = 2, 3).
The response coefficients, κ n , depend on such features of the colliding system as masses of the projectiles, centrality class, or the collision energy, but are to a good approximation [40,41] independent of the eccentricities, hence linearity follows. For higher rank n, as well as for collisions with few participants, non-linear effects [40,42] spoil proportionality (4), hence care is needed in its application.
In practical terms, Eq. (4) means that one can form ratios of flow observables where the response coefficient κ n cancels out, for instance for the cumulant coefficients [42][43][44][45] obtained with k 1 and k 2 particles, or the normalized symmetric cumulants [46,47] obtained with k particles where In addition to NSC( 2 {2}, 3 {2}), in following we frequently use the "double eccentricity ratio" used in [44] as a possible probe of the α clusterization in 12 C+Au collisions. We begin the presentation of our Glauber model results for 16 O + 16 O collisions with the ellipticity and triangularity of the fireball obtained with two-and four-particle cumulants, where specifically    These observables, of course, are not independent of the hydrodynamic response, yet it is worth to have a look at them, as they quantify the shape of the fireball and its fluctuations.
In Figs. 4 and 5 we can see the behavior of the eccentricities. The ellipticity decreases, as expected, with the increasing number of participants, which is the result of the geometry (the fireball is less deformed for the central collisions than for the peripheral collisions). Triangularity, due entirely to fluctuations, at the lower collision energy exhibits a non-monotonic behavior, with maximum around N W = 12.
Passing from the two-to four-particle cumulants reduces the eccentricities, as expected from the consider- ations of fluctuations [43], which increase n {2} compared to n {k}, k = 4, 6, . . . , which are approximately equal [48]. This reduction effect can be noted from Fig. 6 which displays the ratio of eccentricities from four-and two-particle cumulants. In Fig. 7 we investigate the effects of nuclear correlations present in the used distributions from [25] for 2 {2} (solid line), comparing it to the case where the correlations are removed by the mixing technique described in Sec. II (dashed line). We note that the difference is small, at the level of a few percent, with the presence of correlations raising the value of 2 {2} for peripheral collisions, and decreasing it for central collisions. Similar size effects appear for other eccentricity coefficients. In Fig. 8 we present an analogous study of the double eccentricity ratio. We note that the effect of correlations shows up only for the most central events (c < 10%) and reaches of a relative size of about 5% at c = 1%. An analogous study of the normalized symmetric cumulant shown in Fig. 9 leads to a similar conclusion. We note that NSC( 2 {2}, 3 {2}) remains negative for all central- ities (values of N W ) and exhibits a non-monotonic behavior, both at low and high collision energies, as can be inferred from Fig. 10.
The studies of p+ 16 O reactions at the planned √ s N N = 10 TeV collision energy [14] correspond to interactions in air showers at the proton LAB energy of 50 PeV. In Fig. 11 we show our predictions for the double eccentricity ratio and the normalized symmetric cumulant. We note that the behavior is qualitatively similar to the case of 16 O + 16 O from Fig. 10.
In Fig. 12 we show a quantity relevant to cosmic airshower considerations, namely, the p+ 16 O production cross section, plotted against the NN inelastic cross section, which depends on the collision energy. The solid line in the figure presents results within the collision energy range √ s N N = 5 GeV -57 TeV implemented in GLISSANDO 3, whereas the dashed line is an extrapolation made according to the fit formula In Ref. [52], the Pierre Auger Collaboration obtained for √ s N N = 57 TeV the value σ prod p+air = 505 +28 −36 mb with the corresponding inelastic protonproton cross section σ inel p+p = 92 +9 −11 mb. From our fit with Eq. (10) to the p + 16 O GLISSANDO 3 simulations, we get σ prod p+air (92 mb) 507 mb, in a good agreement with [52].  FIG. 14. Same as in Fig. 13 but for the normalized symmetric cumulants (correlated 16 O distributions only).

VI. 16 O REACTIONS WITH HEAVY TARGETS
In our previous papers [44,53,54] we have argued that the heavy-light collisions may reveal the cluster correlations in the light projectile. The best case here is probably the 12 C nucleus, which is believed to have a significant triangular α-cluster component in the ground-state wave function.
We note from Fig. 13 that the double eccentricity ratio is sensitive to the nuclear correlations for the most central collisions, with effects of a relative size of about 10% at c = 1%. For the uniform case, the curves visibly flatten for the most central collisions, whereas with correlations present they continue growing. The behavior is similar at both collision energies.
For the symmetric cumulants shown in Fig. 14 there is only some moderate difference between the correlated and uniform 16 O distributions, hence we present only the correlated case. We note a characteristic non-monotonic behavior with a minimum at low N w and a maximum at intermediate N w .

VII. 16 O + 16 O COLLISIONS WITH WOUNDED QUARKS
In this section we investigate the possible role of nucleon substructure in flow signature of 16 O + 16 O collisions. We compare the predictions of the wounded nucleon model and the wounded parton (wounded quark) model [55][56][57][58] with three constituents, which has turned out successful phenomenologically in explaining the RHIC and LHC data [59][60][61][62][63][64][65][66][67]. The wounded parton picture is implemented in GLISSANDO 3 [18] by placing three partons around the center of each nucleon with an appropriate exponential distribution. The parton-parton inelasticity profile is adjusted in such a way that the resulting NN inelasticity profile generated in p-p collisions matches the phenomenological form discussed in Sec. III.
The comparison of the two models for 16 O + 16 O is shown in Figs. 15 and 16 for the double flow ratio and the normalized symmetric cumulant, respectively. We note that the differences for the double ratio are small, at the level of a few percent, hence the model predictions for these observables are robust with respect to the inclusion of the partonic substructure. For the case of the symmetric cumulant the differences are more visible, with a more prominent minimum occurring for the partonic case.  close collision energy of √ s N N = 2.76 TeV, where also the data have been collected [49,51].
Our results for the double eccentricity ratio are shown in Fig. 17. We note the same pattern in the dependence on N w in all the three reactions, with the minima shifted to the left and upwards with the decreasing projectiles' mass. This feature simply reflects the increasing value of N w with the mass.
For the case of the normalized symmetric cumulant presented in Fig. 18, a similar pattern is observed for all considered collision systems. We note a hallmark nonmonotonicity, and negative values of the function at all values of N w . The cases for Xe-Xe and Pb-Pb collisions agree reasonably well with the data for the symmetric cumulants corresponding to the harmonic flow.

IX. CONCLUSIONS
We have provided a comprehensive Glauber Monte Carlo analysis of ultra-relativistic reactions with 16 O nuclei, including 16 O+ 16 O, p+ 16 O, and 16 O collisions on heavy targets. Although our study is limited to the properties of the initial condition, relying on eccentricities evaluated in the model, it bares significance for experimental studies using harmonic flow, since we apply specific measures approximately independent of the hydrodynamic or transport response of the system. We have also studied the case of the wounded quark model, which leads to similar predictions as the wounded nucleon model.
In our analysis we have used correlated nuclear distributions for 16 Fig. 17 but for the normalized symmetric cumulant. The experimental data come from [50,51].
where the nuclear correlations lead to effects for most central collisions at the level of 10%.
Our basic conclusions are that within the applied collective framework no major qualitative differences should be expected from comparisons of flow characteristics in 16 O+ 16 O collisions to the case of the earlierstudied heavy-ion collisions, such as Xe+Xe or Pb+Pb. If confirmed experimentally, it would hint to a similar collectivity-based mechanism of the fireball evolution across these systems, from small to large. Some opportunities would also be offered by 16 O collisions on a heavy target, where the internal correlation structure in the light nucleus is expected to be of relevance to the harmonic flow characteristics.