# Simulation of neutron background for a dark matter search experiment at JUSL

S. Banik11footnotetext: Corresponding author V. K. S. Kashyap  S. Ghosh  S.
Dutta  B. Mohanty  Meghna K. K
P. Bhattacharjee  S. Saha

###### Abstract

Dark matter search experiments demand low to ultralow radiation background to
operate. It is very important to understand the nature of the radiation
background including knowledge about the sources contributing to it.
Sometimes, evaluation of the background becomes very specific to the site
chosen for the experiment, and also to the experimental configuration. A dark
matter search experiment is proposed to be set up at the Jaduguda Underground
Science Laboratory (JUSL) in India. The laboratory will be located inside an
existing mine with 555 m of vertical rock overburden. Neutrons produced from
$(\alpha,n)$ reactions, spontaneous fission of natural radioactive impurities
in the rocks, and also from cosmic muon induced reactions are considered as
the main background which can affect the sensitivity and outcome of the
experiment. In this work, simulations based on GEANT4 are done to understand
both the radiogenic neutron background caused by natural radioactivity of the
surrounding rock and the cosmogenic neutron background due to interaction of
the deeply penetrating cosmic muons with the rock material. The muon flux in
the cavern is obtained to be $4.49(\pm 0.25)\times 10^{-7}$
cm${}^{-2}\thinspace$s-1 and the fluxes of radiogenic and cosmogenic neutrons
above an energy threshold of 1 MeV in the cavern are obtained to be $5.75(\pm
0.58)\times 10^{-6}\thinspace$cm${}^{-2}\thinspace$s-1 and $7.25(\pm
0.65)\times 10^{-9}\thinspace$cm${}^{-2}\thinspace$ s-1 respectively. The
values obtained are comparable with estimates and measurements done for DAMA,
WIPP and dark matter experiments at Boulby mine. The effectiveness of
different shielding materials are also investigated to obtain the best
possible neutron background reduction for a dark matter search experiment at
JUSL. We also estimate the sensitivity of a CsI based detector for Weakly
Interacting Massive Particle (WIMP) dark matter search at JUSL considering the
estimated neutron background.

##  1 Introduction

Observations of the cosmic microwave background anisotropies indicate that our
universe contains only about 5% of normal or luminous matter, about 27% dark
matter (DM), and the rest as dark energy which accounts for the accelerated
expansion of the universe [1]. DM, as the name implies, has the characteristic
feature that it rarely interacts with normal matter, if at all, to produce
measurable signal in a detector. The DM particles can interact with (scatter
off) the electrons or nuclei in the detector medium. The resulting recoiling
nuclei or electrons would dump their energy in the detector material producing
scintillation, ionization and also heat pulses (such as phonons in solid) [2,
3] depending on the chosen detector material.

It is well known that there will also be electron recoils generated by gamma
rays, beta particles and other ambient radiation quanta which will produce
undesirable background signals. Neutron is an important background that
affects the outcome of highly sensitive experiments such as direct search for
DM[2], neutrinoless double beta decay (NDBD) [4], neutrino oscillation
experiments [5] etc. Neutrons undergo elastic and inelastic scattering with
nuclei of the active medium, resulting in nuclear recoil. The same kind of
nuclear recoil may result from the interaction of DM particles. Therefore,
neutron background limits the sensitivity of direct dark matter search
experiments.

Dark matter search at INO (DINO) is a proposed experiment in India on direct
search for dark matter, to be set up at a future facility of the India-based
Neutrino Observatory (INO) [6]. The experiment will be designed to look for
the signature of DM candidates through observation of extremely tiny amount of
nuclear recoil or electron recoil events in a suitable crystalline detector
medium. Inorganic scintillators, such as Cesium Iodide (CsI) and Gadolinium
Galium Aluminium Garnet (${\rm Gd}_{3}{\rm Ga}_{3}{\rm Al}_{2}{\rm O}_{12}$)
are possible detector materials to explore for the DINO experiment. Both the
scintillators manifest good light yield of around 50 photons keV-1 [7, 8]. The
first phase of the experiment involves exploration and evaluation of the
background effects caused by ambient radiation and other sources at the
specific site. To begin with, a small underground laboratory (approximately 5
m $\times$ 5 m $\times$ 2.2 m), named as Jaduguda Underground Science
Laboratory (JUSL), is already established at depth of 555 m ($\sim$1600 mwe
vertical rock overburden) in an existing mine of Uranium Corporation of India
Limited (UCIL), located at Jaduguda in India for exploring the feasibility of
the experiment. The main reason for setting up such experiments at deeper
underground sites is to reduce the direct effects of cosmic secondary
particles. However, most abundant and deeply penetrating component of the
cosmic rays being muons, they can still reach the experimental cavern located
underground producing neutrons and hadronic showers. Neutrons are generated by
spallation reactions of the muons with the underground rock and shielding
materials. Neutrons are also generated by hadronic showers. These neutrons are
known as cosmogenic neutrons. While the energy spectrum and flux of these
neutrons will depend mostly on the depth of the site and the rock composition,
the other sources of neutrons, known as the radiogenic neutrons, caused by the
($\alpha,n$) reactions and spontaneous fission of naturally occurring 238U,
232Th and their daughters from the decay chains, will depend on both the
natural radioactivity and the elemental composition of the surrounding rocks.
Quantitative estimates for the two sources of neutron background are
therefore, highly dependent on the site, specifically on the abundance of U/Th
and the rock composition. Although the relative flux due to cosmogenic
neutrons are a few orders of magnitude less than that of the radiogenic
neutrons, cosmogenic neutrons are more energetic and therefore, penetrate
deeper inside the passive shielding materials, such as lead, copper and
plastic, used in dark matter experiments. It is important to compare these
backgrounds for different passive shielding configurations to arrive at the
best possible configuration for their effective reduction.

The present study attempts to estimate the neutron background at the JUSL site
using the GEANT4 toolkit [9], with the sources of radiogenic neutrons
generated from standard compilations [10] and source of cosmogenic neutrons
generated from muons using Gaisser parametrization [11]. Propagation of these
neutrons inside the shielding materials are also studied to conclude about the
best possible shielding configurations to setup a DM search experiment at the
site. In section 2, the various sources of neutrons in underground
laboratories are discussed. In section 3, the estimation of radiogenic,
cosmogenic and total neutron flux at the site is explained in detail. In
section 4, various shielding combinations to reduce the neutron flux for a
typical direct dark matter search setup assuming a CsI detector are discussed.
The estimated sensitivity of a CsI detector based experiment at JUSL for
direct dark matter search is presented in section 5. We then conclude with the
consolidation of results and how they compare with other experimental
locations around the world.

##  2 Sources of neutron in underground laboratories

Neutrons in underground laboratories are produced in reactions either by
natural radioactivity or by cosmic rays muons.

##### Radiogenic neutrons:

The radiogenic neutrons are mainly produced in $(\alpha,n)$ reactions caused
by the $\alpha$ particles from the U/Th traces present in surrounding rocks
and the detector materials. The neutrons produced from spontaneous fission of
U or Th and neutron induced fission also contribute to radiogenic neutron
flux. The neutron induced fission, being a secondary process, will have
negligible contribution to the radiogenic neutron flux. The expected energy
spectra of the radiogenic neutrons extend up to about 14 MeV. The radiogenic
neutrons from the rock due to ($\alpha,n$) reactions and due to spontaneous
fission of 238U are considered in the present work.

##### Cosmogenic neutrons:

Cosmogenic neutrons are produced in cosmic ray muon induced interactions
inside rock or shielding materials. Muons produce neutrons via following
interactions:

  1. 1.

interaction with nuclei producing nuclear disintegration.

  2. 2.

muon capture by nucleus followed by neutron emission.

  3. 3.

neutron production by hadrons from muon generated showers.

  4. 4.

neutron production by gammas from muon generated electromagnetic showers.

Cosmogenic neutrons have harder energy spectrum with energies upto few GeVs.
These neutrons can reach detector after propagating large distance, away from
the parent muon track. Because of their energy, it is very hard to shield
them. Some of these neutrons can be vetoed if the associated muons give hit in
the muon veto. The neutrons which are produced inside rock cannot be vetoed.

##  3 Estimation of neutron flux using GEANT4

Simulations have been performed using GEANT4 [9] version 10.02. The reference
physics list “Shielding”, which is recommended for shielding applications at
high energies, is used for describing the physics processes [12]. The high
energy part of this physics list is taken from FTFP_BERT physics list. It uses
high precision neutron model. It was originally developed for neutron
penetration studies and ion - ion collision. It has been used for simulation
of high energy calorimetry and underground or low background experiments.
Secondary particle production threshold cuts are set to 0.7 mm for gammas and
$e^{+}/e^{-}$.

###  3.1 Simulation of radiogenic neutron flux

For the calculation of neutron yield from $(\alpha,n)$ reactions, we have used
the composition of rock surrounding JUSL obtained by wet-chemical, radiometric
and Inductively Coupled Plasma - Optical Emission Spectrometry (ICP-OES)
analyses [13] (given in Table 1). The rock sample was collected by core
drilling at 555 m depth. ICP-OES and wet-chemical analyses give the elemental
composition in the rock whereas the radiometric analysis only gives an
estimate of the K, U and Th content of the rock. Due to the variation among
rock samples and the results obtained from different methods, the uncertainty
in the elemental composition is around 10-20%. The rock, having a density
$\sim 2.89\pm 0.09$ g cm-3, contains 8 ppm of U and 16 ppm of Th [13].

Table 1: The composition of Jaduguda rock as obtained by wet-chemical, radiometric and ICP-OES analyses [13]. Element | Conc (%) | Element | Conc (%)
---|---|---|---
U | 0.0008 | Na | 1.2
Th | 0.0016 | K | 2.2
40K | 0.00034 | Ti | 0.34
Si | 31.0 | P | 0.079
O | 47.8 | Mn | 0.023
Al | 9.6 | Mo | 0.002
Fe | 3.8 | H | 0.028
Ca | 1.3 | S | 0.3
Mg | 0.83 | Others | $<1.5$

The neutron yield $Y_{i}(E_{n})$ from ($\alpha,n$) reactions for a thick
target $i$, have been taken from [10]. The total neutron yield from rock is
calculated by adding the individual neutron yield from the elements weighted
by their mass ratio in the rock. The rock composition information given in
Table 1 is used for this calculation222An intuitive program developed by the
authors of Ref. [10] available at http://neutronyield.usd.edu/ has been used
to obtain the neutron yield by providing the Jaduguda rock composition as
input.. The neutron yields from ($\alpha,n$) reactions in the rock have been
obtained as $6.77\pm 1.12$ yr-1g-1 from 238U and $5.33\pm 0.90$ yr-1g-1 from
232Th.

The specific neutron activity due ($\alpha,n$) reaction in 1 g of rock
normalized to the U and Th content is given in Figure 1. It can be noticed
that the neutrons from ($\alpha,n$) reactions have energies less than 12 MeV.

The natural abundance of 238U is 99.27% with a spontaneous fission half-life
of $8.2\times 10^{15}$ years [14, 15]. The spontaneous fission half-life of
232Th is around 6 orders of magnitude higher and can be neglected [14]. To
simulate neutrons due to spontaneous fission, the Watt function has been used.
The Watt function was initially used to explain fission due to thermal
neutrons in 235U [16]. But it holds good for spontaneous fission of other
heavy nuclei as well. The Watt function is given as

$W(a,b,E^{\prime})=a\sqrt{\frac{4a}{\pi
b}}\exp\left(-\frac{b}{4a}-aE^{\prime}\right)\sinh(\sqrt{bE^{\prime}}),$ (3.1)

where $a=1.54245$ MeV-1 and $b=6.81057$ MeV-1 are constants (for 238U) and
$E^{\prime}$ is the secondary neutron energy [17]. The average neutron
multiplicity ($\bar{\nu}$) for spontaneous fission of 238U is 2.01 [17].
Neutron yield due to spontaneous fission of 238U in the rock is estimated to
be $3.43\pm 0.55$ yr-1g-1.

Figure 1: Expected specific neutron activity of Jaduguda rock due to
($\alpha,n$) reactions and spontaneous fission of 238U.

####  3.1.1 Transmission of radiogenic neutrons through rock

The JUSL has a vertical rock overburden of 555 m. Radiogenic neutrons that are
produced a few meters away from the rock-laboratory boundary, do not
contribute to the neutron flux in the lab. A study is done to find the rock
element thickness to be considered for the calculation. This will help in
reducing statistical fluctuations and efficient calculation of neutron flux.
The rock composition given in Table 1 is used for defining rock in GEANT4. The
radiogenic neutron energy distribution for simulation is generated using the
specific neutron activity given in Figure 1. To calculate the neutron
transmission probability, $10^{5}$ neutrons are randomly generated on the
surface of a rock slab in an area 0.5 m $\times$ 0.5 m with momentum along
$-Z$ direction. Thickness of the rock is varied and the length and breadth are
fixed to 1 m.

Figure 2: The rock slab model used in GEANT4 to calculate the transmission
probability of neutrons.

As neutrons propagate through the rock, they lose energy and get absorbed or
scatter off. Neutrons coming out of the other side of the rock are recorded.
The transmission probability, which is the ratio of the number of output
neutrons to the number of input neutrons is calculated. It is shown as a
function of rock thickness in Figure 3.

Figure 3: Radiogenic neutron transmission probability as a function of rock
thickness.

Only about 0.5% of total neutrons have been transmitted through rock of
thickness 1.0 m. Therefore, for our simulation estimates, we consider 1 m of
rock thickness surrounding the laboratory cavern as active material
contributing to the radiogenic neutron flux.

####  3.1.2 Flux of radiogenic neutrons in JUSL

We consider the cavern to be a hollow cube of outer side 8 m and inner side 4
m. In the cavern, the volume made up of rock is called the ‘Outer cavern’ and
the hollow volume is called the ‘Inner cavern’. To estimate the neutron flux
in the laboratory, contributions from ($\alpha,n$) reactions and spontaneous
fission are taken into account. In the rock, neutron events with the total
energy distribution given in Figure 1 and isotropic angular distribution, are
generated in the 1 m thick grey region of the Inner cavern as shown in Figure
4(a). The energy distribution of 107 neutron events is shown in Figure 4(b).
Some neutrons propagate through rock and some of them reach the experimental
setup. These neutrons while propagating through rock can further produce
neutrons.

(a) (b)

Figure 4: (a) The side view schematic of the cavern as implemented in GEANT4
to calculate the radiogenic neutron flux. (b) Energy distribution of
radiogenic neutrons generated for $10^{7}$ events using the specific neutron
activity shown in Figure 1.

Neutrons reaching the Inner cavern are recorded. The flux of radiogenic
neutrons reaching the laboratory is obtained as 1.12($\pm 0.11)\times
10^{-5}\thinspace$cm${}^{-2}\thinspace$s-1 above 100 keV (mean energy of 1.34
MeV) and 5.75($\pm 0.58)\times 10^{-6}\thinspace$cm${}^{-2}\thinspace$s-1
above 1 MeV (mean energy of 2.18 MeV). The uncertainty shown in the
parentheses include both statistical and systematic uncertainties. The same
convention is followed throughout the paper unless otherwise specified. The
systematic uncertainty arises due to the variation of rock density and is
around 10 %. The flux of radiogenic neutrons reaching the laboratory as a
function of energy is shown in Figure 5.

Figure 5: Flux of radiogenic neutrons reaching laboratory shown as a function
of energy.

###  3.2 Simulation of cosmogenic neutron flux

Cosmic ray muons are created in the upper atmosphere from the interaction of
primary cosmic rays with the atmospheric nuclei. They are the most abundant
charged particles at sea level with an average energy of about 4 GeV and
intensity $\sim$1 cm${}^{-2}\thinspace$min-1 [18]. Since they are minimum
ionizing particles, they can penetrate the rock overburden and reach the
cavern. Cosmic ray muon events can be tagged using active veto shielding, but
stopping them completely is not possible. Even though the flux of cosmic muons
is small in underground laboratories, high energy muons can penetrate the rock
and shielding materials and generate hadrons. Low energy muons get stopped in
the rock.

####  3.2.1 Cosmic muon event generation

To generate realistic neutron flux, the muon flux at the experimental site
needs to be determined. The minimum energy of muons required to reach the
cavern and their maximum lateral displacement need to be calculated.

####  3.2.2 Muon lateral displacement and maximum distance

To reach the cavern from the surface, cosmic muons have to be quite energetic.
Moreover, muons interact with the earth/rock and undergo scattering. Their
initial direction of propagation is altered and their expected position at a
depth is displaced. By calculating the maximum distance traversed by muons of
given energy, we can estimate the minimum energy required by muons to reach
the cavern. Muons with different fixed energies were made to pass through a
cube of rock of side 6 km along the $-Z$ direction. The maximum distance
traversed and lateral displacement were calculated. The results are shown in
Figure 6(a-d). From Figures 6(a) and 6(b), it can be seen that the average
lateral displacement saturates to $\sim 2.3$ m. But the maximum lateral
displacement can be as high as 30 m. From Figures 6(c) and 6(d) we see that
the minimum energy of muon that can reach the cavern (depth 555 m) is around
300 GeV. Therefore, we do not consider cosmic muons of energy less than 300
GeV in our simulation. The lateral displacement helps in making the simulation
faster and is discussed in the next section.

Figure 6: (a) Lateral displacement distribution of muons from their initial
direction of propagation in the rock. (b) Average lateral displacement as a
function of muon energy. (c) Distribution of distance traversed in the rock by
muons of different incident energies. (d) Maximum and average distance
traversed by muons as a function of energy.

####  3.2.3 Calculation of muon flux at the cavern

The latitude, longitude and elevation information of Jaduguda area has been
obtained using Google Earth Pro [19] and the topological profile is shown in
Figure 7. The density of rock is taken to be 2.89 g cm-3 [13].

Figure 7: Elevation map of the area around JUSL (12.2 km $\times$ 5.45 km).
Figure 8: Schematic describing the methodology of calculating muon flux at the
cavern.

The energy and angular distribution of muons on the ground surface are
generated using Gaisser’s parameterization [11]

$\frac{\mathrm{d}^{2}N_{\mu}}{\mathrm{d}E_{\mu}\mathrm{d}\Omega}\approx\frac{0.14E_{\mu}^{-2,7}}{\hbox{cm}^{2}\thinspace\hbox{s}\thinspace\hbox{sr}\thinspace\hbox{GeV}}\times\left[\dfrac{1}{1+\dfrac{1.1E_{\mu}\cos\theta_{z}}{\epsilon_{\pi}}}+\dfrac{\eta}{1+\dfrac{1.1E_{\mu}\cos\theta_{z}}{\epsilon_{K}}}\right],$
(3.2)

where $\theta_{z}$ is the zenith angle, $E_{\mu}$ is the energy of the muon,
$N_{\mu}$ is the number of muons, $\Omega$ is the solid angle,
$\epsilon_{\pi}$ = 115 GeV, $\epsilon_{K}$ = 850 GeV and $\eta$ = 0.054.

This parametrization is valid when the decay of muon is negligible and
curvature of the earth can be neglected (i.e. for zenith angle
$\theta_{z}<70^{\circ}$). The azimuthal angle $\phi$ distribution is generated
uniformly. Muons are propagated to the experimental hall/cavern at 555 m below
the surface. As muon propagates through the earth/rock, it scatters and loses
energy.

Muon events are generated in the range 330 GeV to 15 TeV. Both $\mu^{+}$ and
$\mu^{-}$ events are generated. Since muons impinging on the surface at
$\theta_{z}<70^{\circ}$ are considered, they are generated over a circular
area of radius $\sim$1.5 km with the center at cavern location as shown in
Figure 8. Altitude information ($z$) at each point $(x,y)$ is obtained from
the topological profile (Figure 7). The altitude dependence of muon flux is
ignored. Since a huge number of events333The number of muons that are expected
to fall over a circular area of radius 1.5 km in a day is $\sim 20\times
10^{9}$. are required to be generated on the surface to have a reasonable
number of muons reaching the cavern, the simulation becomes time consuming and
computationally intensive. In order to avoid laborious simulation, we can use
the lateral displacement mentioned in the previous section to save time. From
Figure 6(a), as the maximum lateral displacement is about 30 m, muons with
incident direction within a cone having an axis as the line connecting the
point of incidence to the center of the cavern and radius 30 m (cone opening
angle of $\alpha\sim 3.1^{\circ}$) are simulated (See Figure 8). The other
muons are not simulated and are counted as incident but not reaching the
cavern. The energy distribution of muons before and after reaching the cavern
surface by passing through the earth/rock is shown in Figure 9(a). It is
observed that the flux of muons on the top surface of the outer cavern is
$4.49(\pm 0.25)\times 10^{-7}\thinspace$cm${}^{-2}\thinspace$s-1. The
information of muons such as position, momentum etc. are recorded.

Since the energy of neutrons produced from muon interaction gets attenuated in
the rock, similar to the case of radiogenic neutrons, only a rock element with
finite size will contribute to the flux of cosmogenic neutron in the
laboratory. A study similar to that described in the section 3.1.1 is done for
finding the rock thickness to be considered for simulating the muon induced
neutron flux in the experimental hall.

The rock geometry given in Figure 2 is used for the simulation. The energy
distribution of muons at the surface of the outer cavern shown in Figure 9
(blue histogram) is used. The average energy of muons is around 200 GeV. Muons
are propagated from random positions on a plane of dimension (0.5 m $\times$
0.5 m) (Figure 2(a)) through the rock in the $-Z$ direction. The muon
interactions with rock generate neutrons. The neutrons coming out on the other
side of the rock are recorded. The simulation is repeated for different rock
thicknesses ($t$ = 10 cm, 25 cm, 50 cm, 75 cm, 100 cm, 150 cm).

Figure 9: (a) Energy distribution of muons at the surface and after reaching
the cavern. (b) Neutron produced from muon interaction as a function of rock
thickness. (c) Neutron transmission probability as a function of rock
thickness. (d) The spectra of neutron energy at production
($N_{\mathrm{prod}}$) and after transmission through rock of 150 cm thickness
($N_{\mathrm{out}}$) shown upto 14 GeV.

Number of neutrons produced ($N_{\mathrm{prod}}$) as a function of rock
thickness is given in Figure 9(b). As the rock thickness increases, the
neutron production also increases due to the increase in the probability of
interaction.

####  3.2.4 Transmission of cosmogenic neutrons through rock

Neutron energy spectra at production ($N_{\mathrm{prod}}$) and after
transmission through rock ($N_{\mathrm{out}}$) of thickness 150 cm are shown
in Figure 9(d). The production rate of neutrons for 150 cm rock thickness is
around 0.1 neutron/muon.

The neutron transmission probability which is the ratio of number of neutrons
coming out on the other side of the rock to the number of neutrons produced in
the rock ($R=N_{\mathrm{out}}/N_{\mathrm{prod}}$), as a function of rock
thickness is shown in Figure 9(c). As the rock thickness increases, the
transmission probability decreases. About 7% of total neutrons are transmitted
through rock of thickness 200 cm.

####  3.2.5 Flux of muon induced neutrons in JUSL

To estimate the cosmogenic neutron background at the JUSL, the muon flux
obtained at the surface of the outer cavern is used. The muon events are
generated on the five surfaces of the 2 m thick rock around the cavern as
shown in the Figure 10 (front and back surfaces are not shown in the figure).
Muons are allowed to propagate through the rock and reach the cavern. There
are no muons propagating from the bottom side. While going through the rock,
they generate neutrons and other shower particles like hadrons, gamma and
electrons which then enter the laboratory. Some of the neutrons get absorbed
in the rock itself.

(a)

(b)

Figure 10: (a) Geometry used for simulation. (b) Flux of muon induced neutrons
in the Inner cavern.

The flux of neutrons reaching the cavern is shown in Figure 10(b). It can be
seen that the neutrons produced in rock have energies up to 10s of GeVs. The
muon induced neutron flux in the cavern is found to be $0.93(\pm 0.08)\times
10^{-8}\thinspace$cm${}^{-2}\thinspace$ s-1 with no energy threshold and
$7.25(\pm 0.65)\times 10^{-9}\thinspace$cm${}^{-2}\thinspace$ s-1 above 1 MeV.
The systematic uncertainty is due to the variation of the rock density and the
fact that around 7% of neutrons can still come from distances greater than 2 m
from the cavern boundary.

###  3.3 Total neutron flux in JUSL

Figure 11: Total neutron flux due to radiogenic and cosmogenic sources
expected at the cavern shown as a function of energy.

The total neutron flux from radiogenic origin and muon interactions reaching
the laboratory in the energy range 1 MeV - 100 GeV is shown in the Figure 11.
Above 1 MeV energy threshold, the flux of radiogenic neutrons is 5.75$(\pm
0.58)\times$ 10${}^{-6}\thinspace$cm${}^{-2}\thinspace$s-1, whereas the flux
of neutrons produced by muon interaction in rock is 7.25($\pm 0.65)\times$
10${}^{-9}\thinspace$cm$\thinspace{}^{-2}\thinspace$s-1. For energies less
than $\sim$10 MeV, neutrons flux from radiogenic neutrons is around 3 orders
of magnitude greater than the muon induced neutron flux. For energies above 10
MeV only muon induced neutrons contribute to the spectrum. Therefore the total
neutron flux reaching the cavern/laboratory above 1 MeV energy threshold is
found to be 5.76($\pm 0.58)\times$
10${}^{-6}\thinspace$cm$\thinspace{}^{-2}\thinspace$s-1. Our values are
comparable with neutron flux estimations done for dark matter experiments at
Boulby and WIPP salt mines [20, 21].

##  4 Shielding combinations to reduce neutron flux

A typical experimental setup usually consists of several layers of active and
passive shielding. Active shielding can veto muons and associated neutrons.
Passive shielding systems consist of Lead(Pb) or Iron(Fe) for shielding
gammas, hydrocarbons for moderating neutrons and copper for attenuation of
gammas.

For simplicity, a rectangular geometry of shielding ($x~{}=~{}$1 m, $y=$ 1 m
for all the layers) is considered as shown in Figure 12(a)(side view). There
is an air gap of 1 cm between each layer. Thickness of the materials is
varied. Different combinations of Pb and Polypropylene are studied to find the
optimal shielding composition for neutron reduction.

###  4.1 Reduction of radiogenic neutron flux

Radiogenic neutrons reaching the cavern as obtained in Section 3.1.2, Figure 5
are allowed to pass through the shielding. All the radiogenic neutrons are
stopped by 40 cm thick polypropylene shielding and number of neutrons reaching
Pb surface is found to be zero. Radiogenic neutrons generated from rock have
energy upto a few MeVs and can be shielded using an outer polypropylene layer
only.

###  4.2 Reduction of muon induced neutron flux

The flux of neutrons produced from muon interactions in rock (Figure 12(b)) is
used to generate 104 neutrons on a plane of area (0.5 m $\times$ 0.5 m) which
is at 1 cm above the polypropylene layer (PP1) of thickness 40 cm. Neutrons
crossing the boundary of each layer are recorded.

(a) (b)

Figure 12: (a) Rectangular shielding layers used for simulation. The thicknesses of Pb and PP2 are varied. (b) Energy distribution of generated cosmogenic neutrons. Table 2: Different shielding configurations and their effectiveness. Uncertainties shown are statistical only. | Thicknesses of different |
---|---|---
Configuration | shielding layers (cm) | Transmission (%)
|  PP1 |  Pb |  PP2 |
CFG-1 |  40 |  – |  – | $52.31\pm 0.72$
CFG-2 |  40 |  30 |  – | $136.3\pm 1.17$
CFG-3 |  40 |  30 |  10 | $32.52\pm 0.57$
CFG-4 |  40 |  30 |  20 | $10.44\pm 0.32$
CFG-5 |  40 |  25 |  20 | $12.36\pm 0.35$

About 47% of the input neutrons are stopped by the 400 mm thick polypropylene.
It is observed that a large number of neutrons are produced in Pb from
interactions initiated by incoming neutrons. Neutron back-scattering is also
large at the Pb boundary. A second polypropylene layer (PP2) is needed to
attenuate the neutrons produced in Pb. From Table 2, we see that shielding
combination CFG-4 provides the best neutron reduction. Neutrons, gamma, muons
and electrons are the major backgrounds in a typical experimental setup. Muons
and neutrons can also interact with the detector / shielding materials to
produce more neutrons. We investigate the effect of shielding and calculate
neutron flux at a detector using a simple geometry. The shielding design is
based on the dimensions of shielding materials obtained from CFG-4 shown in
Table 2. The geometry of the rectangular rock element and the experimental
setup is given in Figure 13. The experimental setup consists of a cylindrical
CsI crystal with a radius of 2.2 cm and height 4 cm. Surrounding the crystal
there are cylindrical layers of covering and shielding materials with various
thicknesses: Teflon (0.05 cm), copper (0.6 cm), polypropylene (20 cm, PP1), Pb
(30 cm), polypropylene (40 mm, PP1). There is rock block with thickness of 200
cm surrounding this experimental set up with $\sim 100$ cm of air gap between
them.

Figure 13: Schematic diagram of geometry used for simulation. The crystal,
teflon and copper layers are shown together with wave pattern, PP2 is the
black and light gray checkered region, Pb is the light grey region and PP1 is
the black and white checkered region. The dark grey rectangular region is
responsible for the radiogenic neutron background and black rectangular region
is responsible for the cosmogenic neutron flux.

Muons can interact with the shielding materials to produce high energy
neutrons having energies up to few GeVs. Shielding these neutrons are
difficult. The knowledge of the muon flux is therefore very important for
shielding design. The muon flux reaching the outer boundary of each layer of
the experimental setup obtained from this simulation for 10 days of data are
shown in the Table 3. During passage of cosmic muons through the shield
materials, almost no reduction of muon flux can be seen.

Table 3: The flux of muons at the top surface of different layers. Material | Flux $(\mathrm{cm^{-2}\thinspace s^{-1}}$)
---|---
PP1 | 4.45($\pm$0.24) $\times$ 10-7
Pb | 4.52($\pm$0.26) $\times$ 10-7
PP2 | 4.67($\pm$0.30) $\times$ 10-7
Cu | 3.82($\pm$1.12) $\times$ 10-7

The muons and muon induced neutrons are tracked through the successive layers
of shielding. Neutrons produced in each layer of shielding is recorded.
Neutrons are produced from both muon initiated interactions and neutron
interactions. The neutron production rate in each layer is shown in Table 4.

Table 4: The production rate of neutrons in different layers. Uncertainties shown are statistical only. Material | Neutron production rate ($\mathrm{cm^{-3}\thinspace s^{-1}}$)
---|---
PP1 | 2.05($\pm$0.04) $\times$ 10-10
Pb | 1.72($\pm$0.01) $\times$ 10-8
PP2 | 6.64($\pm$0.06) $\times$ 10-10

Only a fraction of neutrons produced in each layer reach the inner layer of
detector. Others either get absorbed or scatter off. The neutrons which get
reflected back from a layer (i) can be absorbed by the previous layer (ii) get
transmitted out of the setup or (iii) get reflected back again to the same
layer. These effects are taken into consideration to avoid multiple counting.
Neutrons reaching each layer of the detector include neutrons produced in all
previous layers. For instance neutrons reaching the copper layer include
neutrons produced in rock, polypropylene and lead.

The flux of cosmogenic neutrons estimated at the top surface of different
layers for 20 days of data are shown in Figure 14.

Figure 14: Comparison of neutron energy distribution at different layers of
experimental setup

The flux of neutrons that reach the top surface of each layer of the
experimental setup for 20 days of data is shown in Table 5.

Table 5: The flux of neutrons at the top surface of each layer. Uncertainties shown are statistical only. Material |  $E_{\mathrm{neutron}}^{\mathrm{mean}}$ (MeV) | Flux ($\mathrm{cm^{-2}\thinspace s^{-1}}$)
---|---|---
PP1 | 81 | 8.19($\pm$0.11) $\times$ 10-9
Pb | 280 | 3.04($\pm$0.12) $\times$ 10-9
PP2 | 8 | 1.44($\pm$0.02) $\times$ 10-7
Cu | 19 | 7.44($\pm$3.72) $\times$ 10-8
CsI | 9 | 6.15($\pm$4.35) $\times$ 10-8

The increase in neutron flux at the boundary of PP2 is due to the production
of new neutrons in the Pb layer. The increase in mean energy of neutrons at
the Cu layer can be due to the absorption of lower energy neutrons by PP2. The
mean scattering length of neutrons is smaller in hydrogen compared to other
materials like C, Pb or Fe for neutron energies less than $\sim$10 MeV.
Whereas for higher neutron energies, the mean scattering length increases
compared to other materials [22]. Hence, the higher energy neutrons cannot be
moderated easily using hydrogen-based shielding material.

##  5 WIMP sensitivity estimate based on neutron background

We make estimates for sensitivity by following the method suggested in Ref.
[23], and used in KIMS [24] DM search. Considering a dark matter halo model
with a Maxwellian velocity distribution as described in Ref. [23], the total
WIMP event-rate in recoil energy range between $E_{R_{1}}$ and $E_{R_{2}}$ is
given by [24]

$\begin{split}R(v_{E},v_{\mathrm{esc}})&=\frac{k_{0}}{k_{1}}\int_{E_{R_{1}}}^{E_{R_{2}}}dE_{R}\left\\{c_{1}\frac{R_{0}}{E_{0}r}e^{-c_{2}E_{R}/E_{0}r}-\frac{R_{0}}{E_{0}r}e^{-v^{2}_{\mathrm{esc}}/v_{0}^{2}}\right\\},\\\
R_{0}&=5.47\left(\frac{\mathrm{GeV}/c^{2}}{m_{\chi}}\right)\left(\frac{\mathrm{GeV}/c^{2}}{m_{t}}\right)\left(\frac{\sigma_{0}}{\mathrm{pb}}\right)\left(\frac{\rho_{\chi}}{\mathrm{GeV}/c^{2}/\mathrm{cm}^{3}}\right)\left(\frac{v_{0}}{\mathrm{km/s}}\right),\\\
E_{0}&=\frac{1}{2}m_{\chi}v_{0}^{2},\quad
r=\frac{4m_{\chi}m_{t}}{(m_{\chi}+m_{t})^{2}},\end{split}$ (5.1)

where $m_{\chi}$ is the dark matter mass, $m_{t}$ is the mass of a target
nucleus, $\rho_{\chi}=0.3$ GeV cm-3 is the local dark matter density,
$v_{0}=220$ km s-1 is a Maxwell velocity parameter, $v_{\mathrm{esc}}=650$ km
s-1 is the local galactic escape velocity of WIMP, $k_{0}/k_{1}\approx 1$ and
$c_{1},c_{2}$ are constants which depend on the Earth (target) velocity
$v_{E}$, relative to the dark matter distribution as discussed in Ref. [23].
$\sigma_{0}$ is the WIMP-nucleus ‘zero momentum transfer’ cross-section, and
$R_{0}$ is the total event rate (in kg-1day-1) for $v_{E}$ = 0, and
$v_{\mathrm{esc}}=\infty$. Radiogenic neutrons were completely stopped by the
shielding. Only cosmogenic neutrons were able to reach the crystal. From
simulation, the total number of nuclear recoil events due to neutrons within
the energy range 8-60 keV in the CsI crystal is estimated to be $\sim$6
kg-1year-1 (corresponding to a 90% poisson CL of 12 kg-1year-1). The nuclear
recoil energy scale could be converted into the electron equivalent energy
scale using the quenching factors for CsI crystals reported in Ref. [25]. This
turns out to be $\sim$1.5 to 6.5 keV.

The WIMP-nucleon cross-section can be obtained from the WIMP-nucleus cross-
section using the formula [24]

$\sigma_{W-n}=\sigma_{W-A}\frac{\mu_{n}^{2}}{\mu_{A}^{2}}\frac{C_{n}}{C_{A}}$
(5.2)

where $\sigma_{W-n}$ is the WIMP-nucleon cross-section, $\sigma_{W-A}$ is the
WIMP-nucleus cross-section for nucleus of mass number $A$, $\mu_{n}$ and
$\mu_{A}$ are the reduced masses of WIMP-nucleon and WIMP-nucleus systems
respectively and $C_{A}/C_{n}=A^{2}$ for spin independent interaction [24].
Since the CsI crystal has two different nuclei, the combined limit on the
WIMP-nucleon cross-section ($\sigma$) is given by [24]

$\frac{1}{\sigma}=\frac{1}{\sigma_{\mathrm{Cs}}}+\frac{1}{\sigma_{\mathrm{I}}},$
(5.3)

where $\sigma_{\mathrm{Cs}}$ and $\sigma_{\mathrm{I}}$ are the WIMP-nucleon
cross sections limits for Cs and I nuclei respectively. Neglecting the
background contribution from gamma and $\alpha$ particles, the estimated
sensitivity of a CsI based direct dark matter search experiment at JUSL is
shown in Figure 15. Spin independent interaction has been considered and the
WIMP-nucleon cross-section is shown as a function of WIMP-mass. The blue line
corresponds to a 272 g detector as considered in the simulation, running for 1
year. The red dotted line corresponds to a 200 kg detector running for 1 year
assuming the same level of background. Sensitivity estimates in the present
study only assume background events from neutrons; this leads to an optimistic
estimate of the sensitivity. A more realistic calculation will require
consideration of other backgrounds, which are not calculated in this study,
along with neutrons. The parameters such as the quenching factors etc. have
been assumed to be similar to that of KIMS. A better estimate of the
sensitivity will be obtained if all the detector parameters are measured and
understood specifically for DINO. Nevertheless, setting up a direct dark
matter search experiment at JUSL is feasible.

Figure 15: Sensitivity of a CsI based detector at JUSL. Blue line shows the
sensitivity for a 100 kg day exposure (272 g detector for 1 year) and the red
line shows the sensitivity for 200 kg year exposure (200 kg detector for 1
year).

##  6 Conclusion

An estimation of the neutron flux due to radiogenic and cosmogenic sources for
a dark matter search experiment at JUSL, an underground laboratory in India
with 555 m rock overburden, has been done using GEANT4 framework. The
radiogenic neutrons have energies upto few MeVs whereas muon-induced neutron
energy extends up to few 10s of GeVs. A study has also been performed to find
the optimal shielding combination for effective reduction of these neutron
backgrounds.

The specific neutron activity due to ($\alpha,n$) reactions from the
surrounding rock materials has been obtained as $6.77\pm
1.12\thinspace$yr${}^{-1}\thinspace$g-1 of rock from 238U and $5.33\pm
0.90\thinspace$yr${}^{-1}\thinspace$g-1 of rock from 232Th. The specific
neutron activity due to spontaneous fission of 238U is obtained as $3.43\pm
0.55\thinspace$yr${}^{-1}\thinspace$g-1. The flux of radiogenic neutrons
reaching the inner cavern is obtained as $1.12(\pm 0.11)\times
10^{-5}\thinspace$cm${}^{-2}\thinspace$s-1 above 100 keV energy threshold with
a mean energy of 1.34 MeV and $5.75(\pm 0.58)\times
10^{-6}\thinspace$cm${}^{-2}\thinspace$s-1 above 1 MeV energy threshold with a
mean energy of 2.18 MeV.

Cosmic muon events are generated on the surface using Gaisser’s
parametrization and are made to propagate through the rock material. It is
observed that only about 7% of the generated cosmogenic neutrons pass through
the rock of thickness 2 m. The muon flux at the outer cavern is found to be
nearly $4.49(\pm 0.25)\times 10^{-7}$ cm${}^{-2}\thinspace$s-1 with an average
energy of about 200 GeV. These muons are propagated through rock to obtain the
muon and muon induced neutron flux at the cavern. The muon flux at the
boundary of the inner cavern is obtained to be $4.45(\pm 0.24)\times 10^{-7}$
cm${}^{-2}\thinspace$s-1. Muon induced neutron flux from rock in the inner
cavern is found to be $0.93(\pm 0.08)\times 10^{-8}\thinspace$cm-2s-1 without
any energy threshold and $7.25(\pm 0.65)\times
10^{-9}\thinspace$cm${}^{-2}\thinspace$ s-1 above 1 MeV. Our estimated values
of neutron and muon fluxes are comparable with calculations for dark matter
experiments in Boulby mine [20]. The measured value of muon flux at the WIPP
salt mine which is at a similar depth ($\sim 1580$ m.w.e) is $4.77\times
10^{-7}$ cm${}^{-2}\thinspace$s-1 [21]. Our estimation of muon induced neutron
flux is comparable with their calculation (1.6 $\times$
10${}^{-8}\thinspace$cm${}^{-2}\thinspace$s-1) reported in the same paper.

The total neutron flux reaching the Inner cavern/laboratory from both
radiogenic and muon induced reactions in rock is found to be 5.76($\pm
0.58)\times$ 10${}^{-6}\thinspace$cm$\thinspace{}^{-2}\thinspace$s-1 mostly
dominated by radiogenic neutrons. Neutrons produced from muon and neutron
interaction with the shielding materials also contribute to the neutron flux
at the detector. A high number of neutrons are produced in Pb.

Radiogenic neutrons are easily stopped by the 40 cm thick polypropylene layer.
Cosmogenic neutrons can penetrate the shielding and reach the detector.
Moreover, muons generate neutrons while traversing through the shielding. Muon
and neutron fluxes have been estimated at various layers of shielding. It is
seen that a shielding configuration comprising of Polypropylene (40 cm) + Pb
(30 cm) + Polypropylene (20 cm) from the outside towards the experimental
setup gives the best reduction in neutron background. Since the neutron
production rate is high in Pb, a second layer of polypropylene is needed for
effective shielding of these neutrons. With the best determined shielding
configuration, the flux of muon induced neutron at the detector is found to be
6.15($\pm$4.35) $\times$ 10-8 cm${}^{-2}\thinspace$s-1.

The radiogenic neutron yield from contamination in the detector materials has
not been considered in this calculation. D.M. Mei et al. [10] showed that
there is no $(\alpha,n)$ neutron yield in lead due to a very high Coulomb
barrier.

The sensitivity of a CsI based WIMP dark matter search experiment at the
Jaduguda mine has been estimated. The results show that a direct WIMP dark
matter search experiment is feasible at JUSL.

## Acknowledgements

We would like to thank the Department of Atomic Energy (DAE), Government of
India, for support through the project Research in Basic Sciences (Dark
Matter). BM acknowledges the support of JC Bose National fellowship of the
Department of Science and Technology (DST), Government of India. PB and SS
acknowledge the support under Raja Ramanna fellowship of the Department of
Atomic Energy (DAE), Government of India.

## References

  * [1] Planck Collaboration, Ade, P. A. R., Aghanim, N., Arnaud, M., Ashdown, M., Aumont, J. et al., _Planck 2015 results - xiii. cosmological parameters_ , _A &A_ 594 (2016) A13.
  * [2] L. Baudis, _Direct dark matter detection: The next decade_ , _Physics of the Dark Universe_ 1 (2012) 94 – 108.
  * [3] M. Schumann, _Direct detection of WIMP dark matter: concepts and status_ , _Journal of Physics G: Nuclear and Particle Physics_ 46 (Aug, 2019) 103003.
  * [4] C. Alduino, K. Alfonso, F. T. Avignone, O. Azzolini, G. Bari, F. Bellini et al., _Study of rare nuclear processes with CUORE_ , _International Journal of Modern Physics A_ 33 (2018) 1843002, [https://doi.org/10.1142/S0217751X18430029].
  * [5] Daya Bay Collaboration collaboration, F. P. An, A. B. Balantekin, H. R. Band, M. Bishai, S. Blyth, D. Cao et al., _Evolution of the reactor antineutrino flux and spectrum at daya bay_ , _Phys. Rev. Lett._ 118 (Jun, 2017) 251801.
  * [6] A. Kumar, A. M. Vinod Kumar, A. Jash, A. K. Mohanty, A. Chacko, A. Ajmi et al., _Invited review: Physics potential of the ICAL detector at the India-based Neutrino Observatory (INO)_ , _Pramana_ 88 (Apr, 2017) 79.
  * [7] L. Han, B. W. Miller, H. B. Barber, V. V. Nagarkar and L. R. Furenlid, _A new columnar CsI(Tl) scintillator for iQID detectors_ , in _Proceedings of SPIE–the International Society for Optical Engineering_ , no. 9214, 92140D, 2014.  DOI.
  * [8] M. Tyagi, F. Meng, M. Koschan, S. B. Donnald, H. Rothfuss and C. L. Melcher, _Effect of codoping on scintillation and optical properties of a Ce-doped Gd 3Ga3Al2O12 scintillator_, _Journal of Physics D: Applied Physics_ 46 (oct, 2013) 475302.
  * [9] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce et al., _Geant4 - a simulation toolkit_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 506 (2003) 250 – 303.
  * [10] D.-M. Mei, C. Zhang and A. Hime, _Evaluation of ( $\alpha$,n) induced neutrons as a background for dark matter experiments_, _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 606 (2009) 651 – 660.
  * [11] Particle Data Group collaboration, M. Tanabashi, K. Hagiwara, K. Hikasa, K. Nakamura, Y. Sumino, F. Takahashi et al., _Review of particle physics_ , _Phys. Rev. D_ 98 (Aug, 2018) 030001.
  * [12] http://www.slac.stanford.edu/comp/physics/geant4/slac_physics_lists/shielding/physlistdoc.html.
  * [13] “Jaduguda rock analysis results.” Atomic Minerals Directorate for Exploration and Research, India (Private communication.).
  * [14] G. Audi, O. Bersillon, J. Blachot and A. Wapstra, _The Nubase evaluation of nuclear and decay properties_ , _Nuclear Physics A_ 729 (2003) 3 – 128.
  * [15] K. Thiel and W. Herr, _The 238U spontaneous fission decay constant re-determined by fission tracks_, _Earth and Planetary Science Letters_ 30 (1976) 50 – 56.
  * [16] B. E. Watt, _Energy spectrum of neutrons from thermal fission of ${\mathrm{u}}^{235}$_, _Phys. Rev._ 87 (Sep, 1952) 1037–1041.
  * [17] J. M. Verbeke, C. Hagmann and D. Wright, _Simulation of neutron and gamma ray emission from fission and photofission_ , Tech. Rep. UCRL-AR-228518, Lawrence Livermore National Laboratory, Jan, 2014.
  * [18] Particle Data Group collaboration, M. Tanabashi et al., _Review of Particle Physics_ , _Phys. Rev. D_ 98 (2018) 030001.
  * [19] https://www.google.com/intl/en_in/earth/versions/#earth-pro.
  * [20] M. Carson, J. Davies, E. Daw, R. Hollingworth, V. Kudryavtsev, T. Lawson et al., _Neutron background in large-scale xenon detectors for dark matter searches_ , _Astroparticle Physics_ 21 (2004) 667 – 687.
  * [21] E.-I. Esch, T. Bowles, A. Hime, A. Pichlmaier, R. Reifarth and H. Wollnik, _The cosmic ray muon flux at WIPP_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 538 (2005) 516 – 525.
  * [22] C. Bungau, B. Camanzi, J. Champer, Y. Chen, D. Cline, R. Luscher et al., _Monte carlo studies of combined shielding and veto techniques for neutron background reduction in underground dark matter experiments based on liquid noble gas targets_ , _Astroparticle Physics_ 23 (2005) 97 – 115.
  * [23] J. Lewin and P. Smith, _Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil_ , _Astroparticle Physics_ 6 (1996) 87 – 112.
  * [24] H. Lee, H. Bhang, J. Choi, I. Hahn, D. He, M. Hwang et al., _First limit on WIMP cross section with low background CsI(Tl) crystal detector_ , _Physics Letters B_ 633 (2006) 201 – 208.
  * [25] H. Park, D. Choi, J. Choi, I. Hahn, M. Hwang, W. Kang et al., _Neutron beam test of CsI crystal for dark matter search_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 491 (2002) 460 – 469.
  * [26] J. Klinger and V. A. Kudryavtsev, _Muon-induced neutrons do not explain the DAMA data_ , _Phys. Rev. Lett._ 114 (Apr, 2015) 151301.


# Contribution to differential $\pi^{0}$ and $\gamma_{\mathrm{dir}}$
modification
in small systems from color fluctuation effects

Dennis V. Perepelitsa<EMAIL_ADDRESS>University of Colorado Boulder
Boulder, CO 80309

###### Abstract

A major complication in the search for jet quenching in proton– or
deuteron–nucleus collision systems is the presence of physical effects which
influence the experimental determination of collision centrality in the
presence of a hard process. For example, in the proton color fluctuation
picture, protons with a large Bjorken-$x$ ($x\gtrsim 0.1$) parton interact
more weakly with the nucleons in the nucleus, leading to a smaller (larger)
than expected yield in large (small) activity events.
A recent measurement by PHENIX compared the yield of neutral pion and direct
photon production in $d$+Au collisions, under the argument that the photon
yields correct for such biases, and the difference between the two species is
thus attributable to final-state effects (i.e., jet quenching). The main
finding suggests a significant degree of jet quenching for hard processes in
small systems.
In this paper, I argue that the particular photon and pion events selected by
PHENIX arise from proton configurations with significantly different
Bjorken-$x$ distributions, and thus are subject to different magnitudes of
modification in the color fluctuation model. Using the results of a previous
global analysis of RHIC and LHC data, I show that potentially all of the pion-
to-photon difference in PHENIX data can be described by a proton color
fluctuation picture at a quantitative level before any additional physics from
final-state effects is required. This finding reconciles the interpretation of
the PHENIX measurement with others at RHIC and LHC, which have found no
observable evidence for jet quenching in small systems.

††preprint: APS/123-QED

##  I Introduction

Measurements of hard processes in relativistic proton- or deuteron-nucleus
($p/d$+A) at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron
Collider (LHC) have a number of scientific purposes [1, 2], including a
precision determination of the parton densities in nuclei (see [3] for a
recent example), constraints on dynamical processes in the initial state of
the cold nucleus (such as those observed at lower energies [4]), and studies
of the general interplay between soft and hard processes in collisions
involving nuclei. Another key interest, given the robust experimental
signatures of Quark Gluon Plasma-like behavior in these “small” systems [5, 6,
7, 8, 9], is the search for evidence of final-state interactions between the
hard-scattered parton and the dilute system. If the parton shower is modified
in the final state as it propagates through the produced system, it may lead
to a number of effects, the most direct of which is the decreased production
of jets or hadrons at fixed transverse momentum ($p_{\mathrm{T}}$), i.e., jet
quenching [10, 11, 12].

By analogy to large, nucleus-nucleus collision systems, one wants to search in
the most “central”, or highest-activity, collisions where the produced system
is largest and longest-lived. Unfortunately, the experimental selection on
centrality in small systems is challenging due to strong auto-correlation
biases and non-trivial physics effects, with magnitudes larger than that
required for a precision constraint on jet quenching effects. For events with
a generic hard process, the selection is sensitive to upward multiplicity
fluctuations (multiplicity vetoes), leading to an apparently enhanced yield in
central (decreased yield in peripheral) collisions compared to that expected
from the estimated number of nucleon–nucleon collisions,
$\left<N_{\mathrm{coll}}\right>$, in the events [13, 14, 15, 16].

As a separate phenomenon, in extreme kinematic regions characterized by
Bjorken-$x$ $\gtrsim 0.1$, the pattern reverses, with a significant depletion
(enhancement) in the high-activity central (low-activity peripheral) events,
as directly observed in measurements by ATLAS [17], CMS [18], STAR [19], and
PHENIX [20] in $p/d$+A events. In a recent measurement by ATLAS [21], the
produced dijet pair was used to estimate the parton-level kinematics in each
event, strongly suggesting that the modifications follow a universal pattern
in the Bjorken-$x$ of the proton, with their magnitude systematically
increasing with $x$. A particular physics interpretation of these observations
in data is the QCD color fluctuation model [22, 23], which proposes that
proton configurations with a large-$x$ parton interact more weakly with
nucleons in the nucleus than configuration-averaged protons, thus leading to a
decreased centrality signal and the particular pattern of $x$-dependent
modifications described above.

To get around these challenges, experimental searches of jet quenching in
$p/d$+A collisions have therefore been performed in centrality-averaged
(minimum-bias) collisions [24, 25, 26], or using techniques such as jet–hadron
correlations [27, 28] and di-jet asymmetries [18, 19] in central events, all
giving null results for a detectable jet quenching-like signature. A recent
PHENIX measurement [29] in 200 GeV $d$+Au collisions used a different
technique and reported the relative yields of neutral pions ($\pi^{0}$’s) to
those of direct photons ($\gamma_{\mathrm{dir}}$’s), under the argument that
both are subject to the same centrality bias effects, and thus any difference
in the degree of modification should be attributed to centrality-dependent
final-state effects (i.e., which would not be felt by the photons). The PHENIX
measurement reports the quantity $R_{\mathrm{dAu,EXP}}$, defined as the double
ratio of the nuclear modification factors for $\pi^{0}$’s and
$\gamma_{\mathrm{dir}}$,
$R_{\mathrm{dAu,\pi^{0}}}/R_{\mathrm{dAu,\gamma_{\mathrm{dir}}}}$, i.e.,
conceptually defining $R_{\mathrm{dAu,\gamma_{\mathrm{dir}}}}$ to be unity for
all centrality selections. The measurement reports $R_{\mathrm{dAu,EXP}}$ in
the most central 0–5% of $d$+Au events as $0.77\pm 0.03(\mathrm{stat})\pm
0.13(\mathrm{syst})$, where the latter is dominated by a global uncertainty
due to the normalization of the $\gamma_{\mathrm{dir}}$ measurement. If
interpreted as a jet quenching effect, the data are challenging to understand
given the extensive null results described above. For example, the PHENIX data
would require a $\pi^{0}$ $p_{\mathrm{T}}$ spectrum shift of $\delta
p_{\mathrm{T}}\approx-0.21$ GeV in the most central events at RHIC [30],
corresponding to a fractional energy loss of $\delta
p_{\mathrm{T}}/p_{\mathrm{T}}\approx 1$–$2\%$ in $d$+Au events, whereas the
measurement by ATLAS in Ref. [28] sets an exclusion limit at 90% confidence
level of $\delta p_{\mathrm{T}}/p_{\mathrm{T}}<1.4$% for charged hadrons in
$p$+Pb collisions at the LHC. However, these two measurements probe partons
with different initial energies and collisions with different energies and
geometries, and thus may not be in explicit tension.

Notably, the $\pi^{0}$ and $\gamma_{\mathrm{dir}}$ were measured in the same
$p_{\mathrm{T}}$ range of $7.5$–$18$ GeV which, since the $\pi^{0}$ carries
only a fraction of the tree-level parton $p_{\mathrm{T}}$, actually
corresponds to different Bjorken-$x$ ranges. Thus, the measurement technique
defined by PHENIX is directly sensitive to the $x$-dependent color fluctuation
effects described above, in addition to any jet quenching effects. In
particular, the strength of the effect changes quickly with $x$ in the region
of the measurement, and is thus only partially “calibrated out” with the
$\gamma_{\mathrm{dir}}$ baseline. In this paper, I argue that the color
fluctuation effect, as calculated using the results of Ref. [23] with no
additional model modifications or re-tuning to new data, can potentially
explain the full centrality dependence of the observable defined by PHENIX.
When accounting for it, the possible magnitude of any jet quenching effect is
significantly smaller, reconciling the interpretation of the PHENIX
measurement with the established constraints on jet quenching by the other
RHIC and LHC experiments.

##  II Method

The Pythia Monte Carlo event generator [31] was used to simulate the
nucleon–nucleon sub-collisions in $d$+Au collisions which produce $\pi^{0}$’s
and $\gamma_{\mathrm{dir}}$’s in the PHENIX kinematic selections [29], and
determine the distribution of Bjorken-$x$ values for these events. The
specific selections are $p_{\mathrm{T}}=7.5$–$18$ GeV and
$\left|\eta\right|<0.35$ for both species. Pythia version 8.307 was used to
generate events at $\sqrt{s}=200$ GeV, with a mixture of proton–proton,
proton–neutron, and neutron–neutron collisions appropriate for $d$+Au. Pythia
was configured with all HardQCD processes and a $\hat{p}_{\mathrm{T,min}}$
(minimum parton–parton transverse momentum exchange) of 5 GeV for the
$\pi^{0}$ case. For the direct photon case, Pythia was configured with all
PromptPhoton processes and photons which are the result of hadron decays were
rejected (i.e., only those radiated by a quark were selected). This
$\hat{p}_{\mathrm{T,min}}$ threshold was checked that it does not introduce a
kinematic bias in the final-state $p_{\mathrm{T}}$ range of interest. One of
the beam hadrons was chosen to represent the proton or neutron in the deuteron
of a $d$+Au collision and its Bjorken-$x$ value, defined in Pythia as the $x$
at which the PDF in the beam is defined [31], was recorded.

Figure 1:  Distribution in Pythia8 of the values of Bjorken-$x$ in the
deuteron in 200 GeV $d$+Au collisions, for events producing $\pi^{0}$’s
(black) or $\gamma_{\mathrm{dir}}$’s (red) in the kinematic selection used by
PHENIX in Ref. [29].

Figure 1 shows the distribution of Bjorken-$x$ values in the deuteron for
events with a $\pi^{0}$ or $\gamma_{\mathrm{dir}}$ within the PHENIX kinematic
selection. For $\gamma_{\mathrm{dir}}$, the distribution is peaked near $x\sim
7.5/100$ = 0.075, where $7.5$ GeV and $100$ GeV are the minimum
$\gamma_{\mathrm{dir}}$ $p_{\mathrm{T}}$ threshold and the beam energy,
respectively, with a small tail of contributions from events with higher $x$.
For $\pi^{0}$, the distribution has a significantly harder tail to large-$x$
values, reflecting the fact that they are produced by the fragmentation of
hard scattered partons and thus carry only a fraction of their energy.
Notably, the average $x$ value is more than 50% larger for $\pi^{0}$\- than
$\gamma_{\mathrm{dir}}$-producing events, with 28% (9%) of the $\pi^{0}$’s
($\gamma_{\mathrm{dir}}$’s) produced in events with $x>0.2$. In principle, the
magnitude of the color fluctuation effect is expected to be different
depending on whether the high-$x$ parton in the deuteron is a quark or a
gluon. However, in these kinematics, the quark/gluon fraction is similar in
the $\pi^{0}$ and $\gamma_{\mathrm{dir}}$ events and thus the difference in
their $x$ distributions is the main distinction.

The distributions of $x$ values in Fig. 1 were then folded with the results of
the color fluctuation model in Ref. [23] to produce a prediction for the
results of the PHENIX measurement. First, the color fluctuation model paper
primarily reports $R_{\mathrm{CP}}$, the ratio of nuclear modification in
central to peripheral events, which must be converted to the
$R_{\mathrm{dAu,EXP}}$ and other observables reported in the recent PHENIX
paper. This is accomplished using the detailed information provided in an
earlier PHENIX publication [13] on the geometric properties (e.g. the average
number of colliding nucleon–nucleon pairs $\left<N_{\mathrm{coll}}\right>$) of
centrality-selected $d$+Au events. Importantly, while the results here are
presented as a function of the average $\left<N_{\mathrm{coll}}\right>$, the
calculations in Ref. [23] on which they are based were determined using the
distribution of the $N_{\mathrm{coll}}$ values in each centrality bin. Second,
the color fluctuation model presents results in exclusive Bjorken-$x$ bins,
which are then combined in a weighted average with the $x$ distributions for
$\pi^{0}$ and $\gamma_{\mathrm{dir}}$ events above. I highlight that the
parameters of this model, determined using previous RHIC and LHC data, have
not been adjusted or updated before the comparison to the PHENIX data of
interest.

Figure 2:  Calculated nuclear modification factor $R_{\mathrm{dAu}}$ within
the color fluctuation model as a function of the average number of colliding
nucleon pairs, $\left<N_{\mathrm{coll}}\right>$. Calculations are shown for
$\pi^{0}$’s (black) and $\gamma_{\mathrm{dir}}$’s (red) matching the PHENIX
kinematic selection. Since this model does not include any overall
modification, the minimum-bias $R_{\mathrm{dAu}}$ is unity by construction.

Notably, the color fluctuation model in Ref. [23] does not itself include any
centrality-independent modification, i.e., $R_{\mathrm{dAu}}$ for minimum-bias
events is unity by construction. However, there are expected to be small
modifications of the minimum-bias rate for both species due to the nuclear
parton densities [3] and the “isospin effect” for direct photon production
from the valence quark region [32, 33]. Given the specific interest in the
centrality dependence of the observable and the small magnitude of the
minimum-bias-averaged effects compared to the large overall normalization
uncertainties in PHENIX data, I do not model these effects explicitly in this
work. If such effects were included, they would result in a global rescaling
of all calculated $R_{\mathrm{dAu}}$ points in the same direction.

Fig. 2 shows the calculated nuclear modification factors, $R_{\mathrm{dAu}}$,
within the color fluctuation model, applied to the $\pi^{0}$ and
$\gamma_{\mathrm{dir}}$ events in the PHENIX kinematics selection, and within
the PHENIX centrality intervals. Since the minimum-bias $R_{\mathrm{dAu}}$ is
unity by construction (see above), the modifications appear as a relative
enhancement/suppression pattern around unity in peripheral/central events.
Importantly, the magnitude of the relative modifications is stronger for
$\pi^{0}$’s than for $\gamma_{\mathrm{dir}}$’s, reflecting the different $x$
distributions in Fig. 1. Indeed, in Ref. [23], the quantity
$\left<\lambda(x)\right>$, which describes the $x$-dependent weakening of the
inelastic nucleon–nucleon interaction strength, falls rapidly from $\approx
0.8$ to $\approx 0.6$ in this Bjorken-$x$ region, resulting in the difference
seen here.

##  III Comparison to Data

Figure 3:  Comparison of the nuclear modification factor $R_{\mathrm{dAu}}$
for direct photons as a function of $\left<N_{\mathrm{coll}}\right>$, showing
the measurement in PHENIX data (black points) and the calculation from the
color fluctuation model (red horizontal lines). The vertical bars and boxes
around the data points indicate the statistical uncertainties and the overall
16.5% normalization uncertainty, respectively.

Fig. 3 compares the calculated $R_{\mathrm{dAu}}$ for $\gamma_{\mathrm{dir}}$
integrated over the $p_{\mathrm{T}}=7.5$–$18$ GeV range within the color
fluctuation model to that measured by PHENIX, where the quantity is styled as
$(Y^{\gamma_{\mathrm{dir}}}_{\mathrm{dAu}}/Y^{\gamma_{\mathrm{dir}}}_{\mathrm{pp}})/N_{\mathrm{coll}}^{\mathrm{Glauber}}$.
For this comparison, the measured minimum-bias $R_{\mathrm{dAu}}$ for
$\gamma_{\mathrm{dir}}$ is essentially unity, and so no overall rescaling of
the model calculation is applied. The model gives a good description of the
modest centrality dependence of the data, with a relative enhancement and
suppression pattern in going from peripheral to central events. As a reminder,
the $\gamma_{\mathrm{dir}}$ selection typically results in events with
$x\approx 0.11$, and thus the impact of the color fluctuation physics on this
observable is modest in the $x$ range relevant for $\gamma_{\mathrm{dir}}$
production.

Figure 4:  Comparison of the ratio of nuclear modification factors between
that for neutral pions and direct photons,
$R_{\mathrm{dAu}}(\pi^{0})/R_{\mathrm{dAu}}(\gamma_{\mathrm{dir}})$, as a
function of $\left<N_{\mathrm{coll}}\right>$, showing the measurement in
PHENIX data (black points) and the calculation from the color fluctuation
model (red lines). The vertical bars and boxes around the data points indicate
the statistical uncertainties and the overall 16.5% normalization uncertainty,
respectively. The dashed red line assumes that $R_{\mathrm{dAu}}$ in minimum-
bias events for both $\pi^{0}$’s and $\gamma_{\mathrm{dir}}$’s is unity. The
solid red line adjusts the calculation by a factor of 0.9, to better match the
actually measured
$R_{\mathrm{dAu}}(\pi^{0})/R_{\mathrm{dAu}}(\gamma_{\mathrm{dir}})$ in
minimum-bias events by PHENIX.

Fig. 4 compares the measured and calculated “double ratio” of the
$R_{\mathrm{dAu}}$ for $\pi^{0}$’s to that for $\gamma_{\mathrm{dir}}$’s, as a
function of $\left<N_{\mathrm{coll}}\right>$, which is the main result of the
PHENIX analysis. The dashed red line shows the result of the color fluctuation
model, under the assumption that the $R_{\mathrm{dAu}}$ in minimum-bias events
is unity for both species. The calculation gives a good, quantitative
description of the data within its significant normalization uncertainty,
including the slow decrease of the double ratio with increasing
$\left<N_{\mathrm{coll}}\right>$.

In the PHENIX analysis, the minimum bias
$R_{\mathrm{dAu}}(\pi^{0})/R_{\mathrm{dAu}}(\gamma_{\mathrm{dir}})$ is
measured to be $0.92\pm 0.02(\mathrm{stat})\pm 0.15(\mathrm{syst})$. Ref. [29]
notes that this value is consistent with unity within the large overall scale
uncertainty, which is dominated by uncertainties in the $\pi^{0}$ and
$\gamma_{\mathrm{dir}}$ $p$+$p$ references. Thus, the data-model comparison
here has a significant degree of freedom in the overall normalization of the
data. The best match between the model and the central values of the data is
achieved by adjusting the model down by a factor $0.90$, shown as a solid red
line in Fig. 4. I note that this factor is close to the measured minimum-bias
$R_{\mathrm{dAu}}(\pi^{0})/R_{\mathrm{dAu}}(\gamma_{\mathrm{dir}})$ value of
$0.92$, and is well within the normalization uncertainty of the data. This may
alternatively be thought of as scaling the data up by approximately two thirds
of its stated global uncertainty. With this normalization, the model now
provides an excellent description of the PHENIX data in detail. While not
explored here, I note that taking into account the uncertainties in the model
parameters themselves [23] may further improve the data-model agreement. Thus,
no other physics effects, such as a significant centrality-dependent jet
quenching, are needed to describe the PHENIX data.

In the future, several experimental avenues may help to separate the
contribution of the color fluctuation effect from searches for the jet
quenching physics of interest. First, one may design the kinematic selections
with an aim to match the Bjorken-$x$ ranges accessed by the hadron/jet and
photon processes, and thus cancel the impact of the color fluctuation effect.
Second, as motivated in Ref. [34], the two effects are expected to have
different sensitivities to changing the projectile, such as in $p$+Au or
3He+Au data previously recorded at RHIC. This has been performed for
$\pi^{0}$’s in Ref. [35], but not yet for $\gamma_{\mathrm{dir}}$. Third, one
could attempt a centrality selection based on measuring spectator neutrons in
a zero-degree calorimeter (e.g. as used in Refs [27, 28]), which may better
isolate the underlying geometry and not be sensitive to physics effects which
manifest as a bias on the produced multiplicity. Fourth, oxygen–oxygen (O+O)
collisions at RHIC and the LHC [36] would allow for the study of jet quenching
in a small system, in which a single weakly-interacting nucleon in the 16O
nucleus will have a much smaller impact on the overall centrality signal.

##  IV Conclusion

This paper examines a recent PHENIX measurement which found a different degree
of centrality-dependent modification for $\pi^{0}$’s and
$\gamma_{\mathrm{dir}}$’s in $d$+Au collisions at RHIC. The measurement
strategy is ostensibly chosen to calibrate out any centrality-dependent biases
with the $\gamma_{\mathrm{dir}}$ measurement. However, I show that the
particular $\pi^{0}$ and $\gamma_{\mathrm{dir}}$ kinematic selections used in
the measurement select events which arise from different distributions of
Bjorken-$x$ values. As such, the PHENIX measurement is directly sensitive to
the physics effects described by the color fluctuation model, in which
centrality-dependent modifications (which do not arise from jet quenching)
systematically increase with $x$. Using the model parameters determined from
RHIC and LHC measurements reported ten years ago, without any post-diction
updating, I show that a straightforward application of the model to the PHENIX
kinematics gives a good description of the experimental data. After accounting
for the possible role of color fluctuation effects, the evidence for any
remaining final-state effects, such as from jet quenching, is significantly
more limited. This finding reconciles the interpretation of the PHENIX data
with that from other measurements at RHIC and LHC, which have set stringent
limits on the possible amount of parton energy loss in small collision
systems.

###### Acknowledgements.

DVP acknowledges John Lajoie, Riccardo Longo, Jamie Nagle, Peter Steinberg,
and Mark Strikman for useful discussions and input on an early draft of this
paper. DVP’s work is supported by the U.S. Department of Energy, grant DE-
FG02-03ER41244.

## References

  * Salgado  _et al._ [2012] C. A. Salgado _et al._ , Proton-Nucleus Collisions at the LHC: Scientific Opportunities and Requirements, J. Phys. G 39, 015010 (2012), arXiv:1105.3919 [hep-ph] .
  * Salgado and Wessels [2016] C. A. Salgado and J. P. Wessels, Proton–Lead Collisions at the CERN LHC, Ann. Rev. Nucl. Part. Sci. 66, 449 (2016).
  * Eskola  _et al._ [2022] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, EPPS21: a global QCD analysis of nuclear PDFs, Eur. Phys. J. C 82, 413 (2022), arXiv:2112.12462 [hep-ph] .
  * Paul  _et al._ [2022] S. J. Paul _et al._ (CLAS), Observation of Azimuth-Dependent Suppression of Hadron Pairs in Electron Scattering off Nuclei, Phys. Rev. Lett. 129, 182501 (2022), arXiv:2207.06682 [nucl-ex] .
  * Nagle and Zajc [2018] J. L. Nagle and W. A. Zajc, Small System Collectivity in Relativistic Hadronic and Nuclear Collisions, Ann. Rev. Nucl. Part. Sci. 68, 211 (2018), arXiv:1801.03477 [nucl-ex] .
  * Weller and Romatschke [2017] R. D. Weller and P. Romatschke, One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at $\sqrt{s}=5.02$ TeV, Phys. Lett. B 774, 351 (2017), arXiv:1701.07145 [nucl-th] .
  * Adam  _et al._ [2017] J. Adam _et al._ (ALICE), Enhanced production of multi-strange hadrons in high-multiplicity proton-proton collisions, Nature Phys. 13, 535 (2017), arXiv:1606.07424 [nucl-ex] .
  * Sirunyan  _et al._ [2018] A. M. Sirunyan _et al._ (CMS), Elliptic flow of charm and strange hadrons in high-multiplicity pPb collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}=$ 8.16 TeV, Phys. Rev. Lett. 121, 082301 (2018), arXiv:1804.09767 [hep-ex] .
  * Aad  _et al._ [2020] G. Aad _et al._ (ATLAS), Transverse momentum and process dependent azimuthal anisotropies in $\sqrt{s_{\mathrm{NN}}}=8.16$ TeV $p$+Pb collisions with the ATLAS detector, Eur. Phys. J. C 80, 73 (2020), arXiv:1910.13978 [nucl-ex] .
  * Tywoniuk [2014] K. Tywoniuk, Is there jet quenching in pPb?, Nucl. Phys. A 926, 85 (2014).
  * Huss  _et al._ [2021] A. Huss, A. Kurkela, A. Mazeliauskas, R. Paatelainen, W. van der Schee, and U. A. Wiedemann, Discovering Partonic Rescattering in Light Nucleus Collisions, Phys. Rev. Lett. 126, 192301 (2021), arXiv:2007.13754 [hep-ph] .
  * Perepelitsa [2020] D. V. Perepelitsa, Connection between soft and hard probes of small collision systems at RHIC and LHC, EPJ Web Conf. 235, 04003 (2020), arXiv:2005.05981 [nucl-ex] .
  * Adare  _et al._ [2014] A. Adare _et al._ (PHENIX), Centrality categorization for $R_{p(d)+A}$ in high-energy collisions, Phys. Rev. C 90, 034902 (2014), arXiv:1310.4793 [nucl-ex] .
  * Adam  _et al._ [2015] J. Adam _et al._ (ALICE), Centrality dependence of particle production in p-Pb collisions at $\sqrt{s_{\rm NN}}$= 5.02 TeV, Phys. Rev. C 91, 064905 (2015), arXiv:1412.6828 [nucl-ex] .
  * Perepelitsa and Steinberg [2014] D. V. Perepelitsa and P. A. Steinberg, Calculation of centrality bias factors in $p$+A collisions based on a positive correlation of hard process yields with underlying event activity, (2014), arXiv:1412.0976 [nucl-ex] .
  * Loizides and Morsch [2017] C. Loizides and A. Morsch, Absence of jet quenching in peripheral nucleus–nucleus collisions, Phys. Lett. B 773, 408 (2017), arXiv:1705.08856 [nucl-ex] .
  * Aad  _et al._ [2015] G. Aad _et al._ (ATLAS), Centrality and rapidity dependence of inclusive jet production in $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV proton-lead collisions with the ATLAS detector, Phys. Lett. B 748, 392 (2015), arXiv:1412.4092 [hep-ex] .
  * Chatrchyan  _et al._ [2014] S. Chatrchyan _et al._ (CMS), Studies of dijet transverse momentum balance and pseudorapidity distributions in pPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, Eur. Phys. J. C 74, 2951 (2014), arXiv:1401.4433 [nucl-ex] .
  * STA [2024] Event activity correlations and jet measurements in p+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 200 GeV, (2024), arXiv:2404.08784 [nucl-ex] .
  * Adare  _et al._ [2016] A. Adare _et al._ (PHENIX), Centrality-dependent modification of jet-production rates in deuteron-gold collisions at $\sqrt{s_{NN}}$=200 GeV, Phys. Rev. Lett. 116, 122301 (2016), arXiv:1509.04657 [nucl-ex] .
  * Aad  _et al._ [2024] G. Aad _et al._ (ATLAS), Measurement of the Centrality Dependence of the Dijet Yield in p+Pb Collisions at sNN=8.16 TeV with the ATLAS Detector, Phys. Rev. Lett. 132, 102301 (2024), arXiv:2309.00033 [nucl-ex] .
  * Alvioli  _et al._ [2016] M. Alvioli, B. A. Cole, L. Frankfurt, D. V. Perepelitsa, and M. Strikman, Evidence for $x$-dependent proton color fluctuations in pA collisions at the CERN Large Hadron Collider, Phys. Rev. C 93, 011902 (2016), arXiv:1409.7381 [hep-ph] .
  * Alvioli  _et al._ [2018] M. Alvioli, L. Frankfurt, D. Perepelitsa, and M. Strikman, Global analysis of color fluctuation effects in proton– and deuteron–nucleus collisions at RHIC and the LHC, Phys. Rev. D 98, 071502 (2018), arXiv:1709.04993 [hep-ph] .
  * Khachatryan  _et al._ [2017] V. Khachatryan _et al._ (CMS), Charged-particle nuclear modification factors in PbPb and pPb collisions at $\sqrt{s_{\mathrm{N}\;\mathrm{N}}}=5.02$ TeV, JHEP 04, 039, arXiv:1611.01664 [nucl-ex] .
  * Acharya  _et al._ [2022a] S. Acharya _et al._ (ALICE), Nuclear modification factor of light neutral-meson spectra up to high transverse momentum in p–Pb collisions at sNN=8.16 TeV, Phys. Lett. B 827, 136943 (2022a), arXiv:2104.03116 [nucl-ex] .
  * Aad  _et al._ [2023a] G. Aad _et al._ (ATLAS), Charged-hadron production in $pp$, $p$+Pb, Pb+Pb, and Xe+Xe collisions at $\sqrt{s_{{}_{\text{NN}}}}=5$ TeV with the ATLAS detector at the LHC, JHEP 07, 074, arXiv:2211.15257 [hep-ex] .
  * Acharya  _et al._ [2018] S. Acharya _et al._ (ALICE), Constraints on jet quenching in p-Pb collisions at $\mathbf{\sqrt{s_{NN}}}$ = 5.02 TeV measured by the event-activity dependence of semi-inclusive hadron-jet distributions, Phys. Lett. B 783, 95 (2018), arXiv:1712.05603 [nucl-ex] .
  * Aad  _et al._ [2023b] G. Aad _et al._ (ATLAS), Strong Constraints on Jet Quenching in Centrality-Dependent p+Pb Collisions at 5.02 TeV from ATLAS, Phys. Rev. Lett. 131, 072301 (2023b), arXiv:2206.01138 [nucl-ex] .
  * Abdulameer  _et al._ [2023] N. J. Abdulameer _et al._ (PHENIX), Disentangling centrality bias and final-state effects in the production of high-$p_{T}$$\pi^{0}$ using direct $\gamma$ in $d$$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV, (2023), arXiv:2303.12899 [nucl-ex] .
  * Drees [2023] A. Drees, PHENIX Overview (2023), 7th International Conference on the Initial Stages in High-Energy Nuclear Collisions.
  * Bierlich  _et al._ [2022] C. Bierlich _et al._ , A comprehensive guide to the physics and usage of PYTHIA 8.3, SciPost Phys. Codeb. 2022, 8 (2022), arXiv:2203.11601 [hep-ph] .
  * Adare  _et al._ [2013] A. Adare _et al._ (PHENIX), Direct photon production in $d+$Au collisions at $\sqrt{s_{NN}}=200$ GeV, Phys. Rev. C 87, 054907 (2013), arXiv:1208.1234 [nucl-ex] .
  * Aaboud  _et al._ [2019] M. Aaboud _et al._ (ATLAS), Measurement of prompt photon production in $\sqrt{s_{\mathrm{NN}}}=8.16$ TeV $p$+Pb collisions with ATLAS, Phys. Lett. B 796, 230 (2019), arXiv:1903.02209 [nucl-ex] .
  * McGlinchey  _et al._ [2016] D. McGlinchey, J. L. Nagle, and D. V. Perepelitsa, Consequences of high-$x$ proton size fluctuations in small collision systems at $\sqrt{s_{NN}}=$ 200GeV, Phys. Rev. C 94, 024915 (2016), arXiv:1603.06607 [nucl-th] .
  * Acharya  _et al._ [2022b] U. A. Acharya _et al._ (PHENIX), Systematic study of nuclear effects in $p$ $+$Al, $p$ $+$Au, $d$ $+$Au, and 3He$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV using $\pi^{0}$ production, Phys. Rev. C 105, 064902 (2022b), arXiv:2111.05756 [nucl-ex] .
  * Brewer  _et al._ [2021] J. Brewer, A. Mazeliauskas, and W. van der Schee, Opportunities of OO and $p$O collisions at the LHC, in  _Opportunities of OO and pO collisions at the LHC_ (2021) arXiv:2103.01939 [hep-ph] .