<table><caption>**Table 36.8:** Selected experiments using low temperature calorimeters. The table shows currently or soon-to-be operating experiments that will search for dark matter or neutrino properties. The dates refer to the start of the program.</caption><thead><tr><th>Sub-K<br>Experiment</th><th>Location</th><th>Detection<br>mode</th><th>Absorber<br>(Total mass)</th><th>Sensor<br>type</th><th># Sensor<br># Crystal</th></tr></thead><tbody><tr><td colspan="6"><strong>WIMPs</strong></td></tr><tr><td>CRESST III<br>(2016)</td><td>LNGS<br>Italy</td><td>Athermal phonon<br>and scint.</td><td>CaWO<sub>4</sub>/Al<sub>2</sub>O<sub>3</sub></td><td>TES</td><td>10</td></tr><tr><td>EDELWEISS †<br>(_SubGeV)</td><td>LSM Modane<br>France</td><td>Thermal phonon<br>and ion.</td><td>Ge</td><td>NTD Ge</td><td>-</td></tr><tr><td>SuperCDMS<br>(2023)</td><td>SNOLAB<br>Canada</td><td>Athermal phonon<br>and ion.</td><td>Ge/Si</td><td>TES</td><td>24</td></tr><tr><td colspan="6"><strong>Neutrino mass</strong></td></tr><tr><td>ECHo [234]<br>(2012)</td><td>Heidelberg<br>Germany</td><td>Thermal<br>phonon</td><td>Au:<sup>163</sup>Ho<br>(0.2µg)</td><td>MMC</td><td>16</td></tr><tr><td>HOLMES [235]<br>(2015)</td><td>Milan<br>Italy</td><td>Thermal<br>phonon</td><td>implanted <sup>163</sup>Ho<br>(18 µg)</td><td>TES</td><td>1000</td></tr><tr><td>NUMECS† [236]<br>(2015)</td><td>LANL<br>USA</td><td>Thermal<br>phonon</td><td>implanted <sup>163</sup>Ho</td><td>TES</td><td>4096</td></tr><tr><td colspan="6"><strong>0νββ decay</strong></td></tr><tr><td>CUORE [223, 237]<br>(2015)</td><td>LNGS<br>Italy</td><td>Thermal<br>phonon</td><td><sup>nat</sup>TeO<sub>2</sub><br>(741 kg)</td><td>NTD Ge</td><td>988</td></tr><tr><td>CUPID [224]<br>(2015)</td><td>LNGS<br>Italy</td><td>Phonon<br>and scint.</td><td>Li<sub>2</sub><sup>100</sup>MoO<sub>4</sub><br>(7450 kg)</td><td>NTD Ge</td><td>1596</td></tr><tr><td>AMoRe-I [227]<br>(2018)</td><td>Y2L<br>South Korea</td><td>Phonon<br>and scint.</td><td>Ca<sup>100</sup>MoO<sub>4</sub>/Li<sub>2</sub><sup>100</sup>MoO<sub>4</sub><br>(6 kg)</td><td>MMC</td><td>13/5</td></tr></tbody></table>

manifest themselves as electron recoils in the detector.

In recent years, multi-ton liquid noble detectors have outclassed LTD-based technologies in searches for heavy (>10 GeV/c²) dark matter owing to their ability to more easily and cheaply scale to large target masses. However, the lower thresholds achieved by LTDs continue to make them the technology of choice for low-mass dark matter searches. New advances have enabled these detectors to reach much lower energy thresholds than previously obtained, albeit sometimes at the cost of being able to detect energy in more than one channel as described above. For example, the use of an electric field to generate Neganov-Trofimov-Luke [238, 239] phonons in proportion to the applied voltage, has enabled the detection of single electron hole pairs in Si detectors with thresholds as low as a few eV (SuperCDMS HVEV). This and similar advances in Ge LTDs (EDELWEISS\_SubGeV) have enabled sensitive searches for dark photons and dark matter that scatters off electrons [240, 241]. Next generation experiments such as SPICE/HERALD aim to further optimize the intrinsic energy resolution of TES detectors, coupled with a strategic choice of target materials (superfluid He and polar crystals) to enable sensitivites to dark matter with masses below an MeV/c² [242]. Current state-of-the-art axion searches use SQUID based quantum amplifiers such as Josephson Parametric Amplifiers along with resonant cavities operating below 100 mK to look for a signal above fluctuations in the thermal noise [243]. Future axion experiments are also working to close the sensitivity gap between particle and wave-like dark matter with the help of LTDs. Broad-band axion searches in the THz range are being proposed, which will make use of TES, SNSPDs or KIDS for single photon detection [244]. Finally, LTD-based dark matter detectors are also actively employed to study coherent neutrino scattering, owing to the fact that the hypothesized signal from dark matter-nucleus scattering is nearly identical to that from neutrino-nucleus scattering, with both inducing nuclear recoils in a similar energy range [245, 246].

## 36.6 Low-radioactivity background techniques

Revised October 2023 by Al. Ianni (INFN, LNGS) and S. Schönert (Munich Tech. U.).

### 36.6.1 Introduction

The study of rare phenomena in fundamental physics, such as proton decay, neutrinoless double beta decay, dark matter,

and MeV-scale neutrino interactions, requires extremely low levels of background radiation. Experiments searching for these rare events record electron recoils or nuclear recoils in the energy scale from a few eV to several MeV. The detector technologies used are multiple from organic and cryogenic liquid scintillators, to bolometers, solide state calorimeters, gaseous detectors, and crystal scintillators. As far as the background contamination is concerned at some extent the application defines the requirements, although the common denominator is that an extreme reduction of all background sources is essential. Leading experiments in rare events search have obtained in the region-of-interest (ROI) a background as low as ~ $10^{-4}$ cts/kg/y/keV. As a first and crucial step, a dedicated radio-purity assay of the detector set-up components has to be carried out. Over the last fifty years, special screening and cleaning techniques have been developed to measure and mitigate ultra-low levels of background. In order to characterize the background sources we refer to Heusser [247] and identify the following five main categories:

* environmental radioactivity in the location where the detector is installed;

* radioimpurities in the detector and shielding;

* radon and its progenies;

* cosmic rays and induced radioactivity;

* neutrons from natural fission, (α,n) reactions, and from cosmic-ray muons interactions.

The energy range affected by these background sources is mainly <10 MeV. All materials contain traces of long-lived primordial radioimpurities, such as $^{238}$U, $^{232}$Th, $^{235}$U ($^{238}$U/$^{235}$U ~ 138), and $^{40}$K ($^{40}$K/$^{K^{nat}}$ ~ $1.17 \times 10^{-4}$). We recall that 1 ppt of $^{238}$U and $^{232}$Th corresponds to 12.36 and 4.06 µBq/kg, respectively; 1 ppb of $^{K^{nat}}$ corresponds to 30.25 µBq/kg. In the Earth's crust the abundance of uranium and thorium is of the order of 1 – 10 ppm which corresponds to about 10 – 100 Bq/kg. Taking into account these contamination levels and the low background requirements, a fundamental background reduction and mitigation is essential to carry out rare phenomena research.

Besides primordial radionuclides other radioactive elements are produced through interactions with matter of secondary cosmic ray particles. Among these so-called cosmogenic radionuclides we recall, in particular, $^{3}$H, $^{14}$C, $^{7}$Be, $^{39}$Ar, $^{42}$Ar, and, in copper, steel and iron often used as shielding materials, $^{57,60}$Co.

<footnote>†No payload size quoted for experiments that are primarily in R&D phase.</footnote>
A third category of background source in our environment con-
sists of anthropogenic radionuclides. These are artificially pro-
duced mainly through nuclear reactions in nuclear power plants,
nuclear fuel reprocessing plants, and nuclear weapons testing. An-
thropogenic background elements of concern for rare phenomena
are $^{85}$Kr, $^{137}$Cs, $^{241}$Am, $^{60}$Co, and $^{90}$Sr. The concentration of
$^{85}$Kr ($T_{1/2}$=10.76 y, $Q_\beta$=687 keV) in air has been slowly increas-
ing since World War Two with a present activity of order 1 Bq/m$^3$.
As a consequence $^{85}$Kr is a crucial background in experiments
making use of nitrogen, xenon, and argon from air.

Once detector components radioactive backgrounds have been
assayed and reduced by a careful selection campaign, which can
last several years, a further step to push the background level be-
yond the screening possibilities is needed to reach the required
sensitivity. For this purpose special active vetos, purification
methods, and offline analyses have been developed. A meticulous
background understanding and mitigation is crucial to explain
any possible signal excess which may be detected. Background
mitigation techniques are based on:

* use of radio-pure materials that absorb ionizing radiation;

* identify radio-pure material for detector construction;

* perform advanced surface and sub-surface cleaning treatment;

* reduce muon flux with underground detector deployment and using active vetos;

* exploit advanced detection and tagging techniques to discriminate signal from background.

In the following we describe radio-assay and background mitig-
ation techniques developed and exploited in the framework of rare
phenomena searches in deep underground laboratories.

36.6.2 Radio-purity assay

The radio-purity assay of detector components is a basic pre-
requisite to be carried out in low counting experiments. Several
techniques are exploited for radio-purity assay. They are comple-
mentary to characterise the radio-purity of materials for shielding
or core detector components. Next generation experiments ra-
dioassay campaign requires a considerable effort and organization
for several years [248–251].

Gamma spectroscopy (GS) via high purity germanium (HPGe)
detectors is a powerful and crucial technique [247]. It is nonde-
structive and thanks to the energy resolution allows to distinguish
various radionuclides elements. Radiation from $^{238}$U and $^{232}$Th
comes with all decay products in the radioactive chains. How-
ever, if secular equilibrium is broken, this crucial information can
be addressed by separating different gamma-ray lines character-
izing elements in the decay chains. In the $^{238}$U chain one can
have three sub-chains out of equilibrium. The first sub-chain can
be assessed through $^{234}$Th direct progeny of $^{238}$U. The second
sub-chain, which originates from $^{226}$Ra, can be probed by $^{214}$Pb
and $^{214}$Bi. A third sub-chain, which starts from $^{210}$Pb, cannot
be efficiently probed by GS, yet is of crucial importance and al-
ternative methods must be used. In the $^{232}$Th chain again one
can probe two sub-chains which can be out of secular equilib-
rium: the first one through $^{228}$Ac from $^{228}$Ra; the second from
$^{228}$Th can be measured through $^{212}$Pb, $^{212}$Bi and $^{208}$Tl. The
technology for HPGe operated in deep underground counting fa-
cilities [252,253] or in shallow laboratories with an efficient active
veto shielding [254] has been boosted to sensitivities of the or-
der of 10-100 µBq/kg by carefully selecting detector components,
electronics and sample handling systems. The HPGe screening
method requires order of ~ 0.1-a few kg of material, and weeks
of acquisition to produce a reliable measurement.

A second crucial technique is based on inductively coupled
plasma mass spectrometry (ICP-MS) [255]. This technique can
probe primordial parent activities at the level of 1 µBq/kg. It is a
destructive method and often needs special sample preparation on
a small quantity of material. The ICP-MS does not measure the
radioactive decay of isotopes but determines their concentration.
At present, it is the most sensitive and rapid screening technique
which allows to select materials at sub-ppt level of impurities. The
drawback is that ICP-MS cannot assess whether the uranium and

thorium chains are out of equilibrium and to reach ultra-high sen-
sitivities one needs to carefully prepare and handle the samples
in a cleanroom environment. ICP-MS screening must be coupled
with other methods to properly assess the radio-purity of mate-
rials in the context of rare events searches. The Glow Discharge
Mass Spectroscopy (GDMS) is a trace element analysis technique
somehow alternative to ICP-MS. An advantage of GDMS is the
possibility to determine the bulk composition of the sample, as-
assuming homogeneity. Sensitivities of the order of 10 ppt can be
achieved.

A third screening technique is based on neutron activation anal-
ysis (NAA) [256]. A sample exposed to a neutrons flux can be ac-
tivated to form radioactive isotopes which can be detected using
HPGe detectors or ICP-MS. Considering the difficulties to irra-
diate samples, this method is less often used. NAA can probe
sensitivities at the level of 0.01 µBq/kg.

In case out of secular equilibrium conditions are measured for
uranium and thorium from GS screening, a rigorous assessment
cannot avoid radon emanation measurements. This matter is dis-
cussed in Section 36.6.3.

Complementarity between radio-purity assay techniques is a
crucial parameter to design detectors for rare phenomena searches.
As it has been pointed out above, requirements can be more strin-
gent than the best sensitivity which can be obtained with the cur-
rent radioassay techniques. In these cases, a number of prototype
detectors have been built and operated to prove the feasibility to
reach ultra-low backgrounds.

36.6.3 Radon and its progeny

Radon is considered to be rare in nature because most of its
isotopes are short-lived. However, $^{222}$Rn ($T_{1/2}$=3.82 d) is of par-
ticular concern in our context. $^{222}$Rn is produced by $^{226}$Ra and
is a radioactive noble gas which can move within active detector
components. $^{222}$Rn daughters are heavy metals which can deposit
on surfaces. Diffusive $^{222}$Rn, supported by $^{226}$Ra, can deposit on
surfaces the long-lived $^{210}$Pb, which is a major concern for low
counting experiments. In addition, due to ~100 keV nuclear re-
coil energies from alpha decays in the $^{226}$Ra sub-chain, eventually
$^{210}$Pb can be implanted into a sub-surface layer of a material ex-
posed to radon. This sub-surface contamination can remain even
after surface cleaning. Surface contamination of $^{210}$Pb is a seri-
ous background for direct dark matter experiments through alpha
decay of $^{210}$Po, which can generate neutrons by (α,n) reactions.
Low energy beta/gamma emissions from $^{210}$Pb are also a concern.
Therefore, radon-free cleanrooms are essential for the assembly of
the detector components. Effective radon abatement systems are
available for this purpose.

In assembling and commissioning rare events experiments, spe-
cial care must be dedicated to the estimate of radon emanation of
the materials and continuous radon monitoring. For this purpose
different methods for radon assay have been developed and ex-
ploited since the beginning of solar neutrino observations. $^{222}$Rn
atoms are collected inside an exhalation chamber for several half-
lives before adsorption and counting. Detection limits of the or-
der of 100 µBq (about 50 $^{222}$Rn atoms) can be obtained with ~50
l stainless steel electro-polished chambers [257]. This limit can
be pushed down to 30 µBq for 1-liter scale chamber. Emanation
of large vessels (cryostats, storage tanks, purification columns)
can be determined by collecting exhaled radon into transportable
charcoal traps [258]. The same method can be used for liquid
samples. In this case instead of evacuating the exhalation cham-
ber into a charcoal trap, He is flushed through a sparger tube for
about 10 times the volume of the liquid used. Sensitivities of the
order of 10 µBq/kg have been reached.

Nitrogen or synthetic air is often used in rare events exper-
iments for purging, stripping, and assembling the experimental
apparatus. These gases might contain radon. In gases $^{222}$Rn can
also be detected using electrostatic collection of $^{218}$Po and $^{214}$Po
[259, 260]. Sensitivities of the order of mBq/m³ can be obtained.
In Borexino three grades of nitrogen purity were used: regular pu-
licity, high purity, and low argon and krypton purity. The regular
purity is obtained from boil-off gas and has radon, measured with
the method reported above, at the level of < 100µBq/m³. For
stripping the purified liquid scintillator a higher purity is needed.
To remove radon from regular purity nitrogen a dedicated ab-
sorber plant has been built. This system can reduce radon by a
factor of 100. Finally, we mention that not only radon is found in
nitrogen. For specific applications the long-lived $^{39}$Ar and $^{85}$Kr
in nitrogen are an important source of background.

36.6.4 Surface backgrounds

Surface contamination of long-lived $^{222}$Rn daughters can be challenging in low-counting experiments. Considering required sensitivities of the next generation experiments, this source of background has to be properly quantified and mitigated. Therefore, exposure to $^{222}$Rn should be monitored and limited to reduce build-up of $^{210}$Pb on surfaces. In addition, $^{222}$Rn exposure could also produce sub-surface contamination as discussed in Section 36.6.3. For one cannot avoid radon contamination in many circumstances during production of detectors components, it is crucial to quantify the effectiveness of cleaning aiming at removing surface contamination of $^{210}$Pb, $^{210}$Bi, and $^{210}$Po. A simple cleaning procedure, which implies degreasing, wiping, and rinsing the material surfaces is not effective in removing these surface contaminants. Studies of cleaning procedures have been carried out exposing stainless steel, copper and other materials to a strong radon source. Etching and electropolishing with subsequental passivation and rinsing have been investigated in great details. Several recipes for etching and electropolishing have been proposed [261]. Electropolishing has been shown to be very effective in reducing $^{210}$Pb, $^{210}$Bi, and $^{210}$Po from both copper and stainless steel by a factor greater than 100. Etching, which is easier to perform than electropolishing, followed by passivation and rinsing with deionized water is effective in reducing $^{210}$Pb and $^{210}$Bi by a factor between 50 and 100. However, it is less effective for $^{210}$Po, which in copper is very poorly reduced. Removing $^{210}$Po from surfaces is crucial, therefore, naturally $^{210}$Po contaminated copper and stainless steel surfaces have been deeply investigated with a high sensitivity (1 mBq/m$^2$) alpha spectrometer [262]. After multi-etching steps (≥ 3), followed by a passivation, a reduction of order 100 has been obtained. On the contrary, static etching (single step) is poorly effective. Electropolishing, or multi-etching are recommended in case copper is in direct contact with the active core of the detector and ultra-high radio-purity is essential or when copper electroforming cannot be used. Electroplating of a thin layer of high radio-pure copper onto the surface of less radio-pure copper has also been investigate to mitigate surface background [262]. This technique is shown to be efficient in reducing surface activity of $^{210}$Po when electroformed copper is used [262].

Besides copper and steel other materials often used, such as polyethylene and teflon, have been investigated to understand how to reduce radon plate-out contaminations [263].

As far as surface contamination is concerned particulate fallout
in cleanroom environments could be of concern. In general, chem-
ical composition of dust reflects local composition of soil and dust.
This is not necessarily true in cleanroom spaces, where the com-
position of dust depends on ongoing activities and handled mate-
rials [264]. The rate of fallout is an important information in the
framework of rare events experiments, where assembling is per-
formed in cleanrooms. The $^{210}$Pb contamination, inferred from
stable lead measured by ICP-MS, due to dust fallout has been
investigated [264]. This contamination has a different origin with
respect to $^{210}$Pb from radon progeny implantation. Ultimately,
one can conclude that radon exposure is more crucial than dust
fallout. Therefore, the best practice for rare events experiments
to face surface background contamination is to perform cleaning
and assembling in a radon-free cleanroom environment.

36.6.5 Mitigation of backgrounds and active background dis-
crimination

In this Section we discuss a selection of different background
mitigation techniques used in experiments to search for rare events
in deep underground laboratories. This includes both the avoid-
ance or reduction of specific radioactive contamination, as well as
active background suppression techniques based on specific event
features and topologies.

* Mitigation of $^{222}$Rn daughters deposition. In dark matter

direct search radiogenic (α,n) reactions due to radioactive
decays are of great concern. The emitted neutron can mimic
a nuclear recoil induced by a dark matter particle interaction.
In the DarkSide-50 and LUX-ZEPLIN (LZ) experiments to
mitigate this background, cleaning of parts and assembling of
the dual-phase Time Projection Chamber have been carried
out in a radon-free cleanroom with $^{222}$Rn activity of the order
of $\lesssim 50$ mBq/m$^3$.

*   **Underground argon.** Liquid argon is an excellent scintillator to search for dark matter interactions due to the high electron-recoils rejection power through pulse shape discrimination. Events from $\beta-\gamma$ background can be rejected at the level of $10^7$ or better with respect to nuclear-recoils. Liquid argon is very effective as active shield and has been used in experiments for neutrinoless double beta decay. However, a major drawback for dark matter search using argon as a target is due to the fact that atmospheric argon contains about 1 Bq/kg of cosmogenic $^{39}$Ar ($T_{1/2}=269$ y, $Q_{\beta}=565$ keV) [265]. For neutrinoless double beta decay search using argon as an active veto, $^{42}$Ar ($T_{1/2}=32.9$ y, $Q_{\beta}=599$ keV) with its short lived progeny $^{42}$K ($T_{1/2}=12$ h, $Q_{\beta}=3525$ keV) is a major source of background. $^{39}$Ar limits significantly the dark matter and $^{42}$Ar the neutrinoless double beta decay sensitivity search. Therefore, a source of argon with reduced $^{39}$Ar and $^{42}$Ar is crucial. Centrifugation or differential thermal diffusion are established methods to separate $^{39}/^{42}$Ar and $^{40}$Ar. However, this is an expensive method for a large fiducial mass. Argon from underground natural gas reservoirs is shown to contain low $^{39}$Ar [266], and it is expected that $^{42}$Ar is similar or even strongly reduced. Therefore, the use of underground argon mitigates the $^{39}/^{42}$Ar backgrounds. DarkSide-50, with a mass of 150 kg of underground argon, has shown that this source of argon contains $^{39}$Ar at a level reduced by a factor of $(1.4 \pm 0.3) \times 10^3$ with respect to atmospheric argon [267]. DarkSide-20k and LEGEND-1000 are planning to exploit underground argon from the Urania facility with a production capacity of 330 kg/day [268]. The former to reduce $^{39}$Ar and the latter $^{42}$Ar.

*   **Electro-formed copper.** High radio-purity copper is often used as shielding and as core detectors components. The radioim-purities in the copper can be a dominant source of exter-nal background. Copper electro-forming is a technique used to reduce this background component in rare events exper-iments. Copper electro-forming is a well known process to obtain ultra-high radio-purity copper. This technique has been used in the framework of the Majorana experiment to search for neutrinoless double beta decay and in ANAIS to search for dark matter annual modulation. Sub-ppt levels in uranium and thorium have been achieved with electro-formed copper [269].

*   **Background suppression using topological event information.** In addition to rigorous selection of high-purity target and shielding materials from external radiation, additional active suppression techniques must usually be employed in low background experiments to achieve the appropriate experimental sensitivities. While signal and background events may be indistinguishable if only their energy deposition is measured, their event features may differ significantly in time and space. Liquid scintillators use the characteristic photon emission time distributions to distinguish between electron- and alpha-like signals [270] and nuclear recoils [271], respectively. High-purity germanium detectors use the time evolution of the induced charges to separate point-like signal candidates for neutrinoless double-beta decay events from background signals induced, for example, by gamma interactions with multiple interactions within a crystal, or from $\beta$ or $\alpha$ events on the n+ or p+ electrodes [272]. The operation of bare high-purity germanium detectors in an instrumented liquid argon shield enabled the GERDA and LEGEND experiments to identify backgrounds with signal-like event topology within the HPGe detectors, but with characteristic energy deposition in the surrounding liquid argon. These synergistic back-
ground suppression techniques enabled for a first time ever a quasi background-free search for neutrinoless double beta decays with GERDA [273].

EXO-200 [190] exploited topological event information distinguishing single-site against multi-site events with a liquid xenon TPC for neutrinoless double beta decay in $^{136}$Xe to reject background and thereby enhancing signal discovery sensitivity. With a gas xenon TPC the NEXT experiment will fully exploit the differences in the spatial ionization patterns of double beta decay and single electron events [274]. The former is characterized by two Bragg peaks at opposite ends of the tracks, the latter on the contrary displays only one peak. The combination of topology information and good energy resolution offer a powerful tool for background rejection. In general, the 3D reconstruction of detected events in a multivariate fit, which accounts for spacial surface and bulk distributions of signal and background together with other properties, such as pulse shape and topological features, is a powerful tool for background mitigation.

*   **Signal selection in direct Dark Matter experiments.** In direct dark matter experiments electron recoil events have to be mitigated with respect to nuclear recoil events. In semiconductor bolometers, operating at a few tens of mK under a bias electric field and used as calorimeters, drifting charges produce a large phonon signal proportional to the ionization, which allows to discriminate electron recoils by the combination of charge and phonon signals [275]. In scintillating bolometers the phonon and light signals are used for the same purpose [276]. Cryogenic scintillators, such as xenon and argon, in time projection chambers offer a strong electron recoils background mitigation through the detection of a primary scintillation signal in liquid and a secondary signal in gas from the drift and extraction of ionization electrons. This background mitigation technique is also being used for neutrinoless double beta decay with $^{136}$Xe. The accurate fiducialization and good rejection of multiple-scattering events allow dark matter optimized experiments to attempt a search for this very rare phenomenon.

*   **Neutron tagging.** Present and next generation experiments need large neutron tagging detectors. In DarkSide-50 a dedicated active veto has been developed to both suppress and measure *in situ* the rate of neutron-induced background events [277]. The detector consisted of a boron-loaded liquid scintillator, which served both as shielding against γ-rays and as a tag for neutrons. Neutrons are thermalised and captured on $^{10}$B. Experimental data has shown a neutron rejection power greater than 99.1% with 5% concentration of TMB in 30 tonnes pseudocumene-based liquid scintillator. The LZ direct dark matter detector with a central TPC of 7 tonnes of liquid xenon makes use of an outer Gd-loaded liquid scintillator neutron tagging veto, which works similarly to the DarkSide-50 veto, replacing boron with gadolinium. The 22 tonnes liquid scintillator is based on linear alkyl benzene (LAB) as solvent. This detector has been designed to operate with a neutron tagging efficiency greater than 95%. The nuclear-recoil background is reduced by a factor of 10. The XENONnT detector is making use of a cylindrical stainless steel tank filled with Gd-loaded water, which surrounds the cryostat with a TPC of 5.9 tonnes of active liquid xenon. The Gd concentration in water is 0.2% in mass. Neutrons leaving the TPC volume will be moderated and captured by the Gd with a probability of 91%. The gamma-rays emitted after the capture are detected from Cherenkov photons, providing a neutron tagging [278].

Neutron tagging techniques are also crucial for coherent neutrino-nucleus elastic scattering detection from pion-decay-at-rest neutrino source at spallation sources [279,280]. Neutrons can be produced in the Pb shielding from CC and NC interactions and propagate into the nuclear recoil sensitive volume, generating a beam-related background. Careful shielding design and implementation of neutron tagging is crucial.

*   **Mitigation of cosmogenic background.** In recent years the required sensitivity to search for dark matter and neutrinoless double beta decay asks for a strong reduction of cosmogenic background. Muons can produce neutrons that can enter the active volume of the detector from the surrounding rock or from external detector components. The yield of these so-called cosmogenic neutrons depend on the muon energies and the material properties of the medium the muon passes through. In addition, muons can produce by spallation radioisotopes inside the detector active volume. Cosmogenic backgrounds are a function of depth and experimental design, and can limit the sensitivity to search for rare events. Most of these radioisotopes are short-lived and their effect can be easily removed by an active veto based on the time correlation with a crossing muon. However, a number of cosmogenic radioisotopes are long-lived and they require an important consideration. Mitigation of these cosmogenic backgrounds produced in-situ deep underground are a major challenge for upcoming and future experiments. Optimizing the detector design and analysis strategies at a given depth equals an effective muon flux reduction. In particular, we mention two discrimination techniques for muon-induced isotopes: 1) $^{11}$C tagging by a three-fold coincidence between the crossing muon, the capture of the ejected neutron from $^{12}$C, and the $^{11}$C decay [281]; 2) similar delayed coincidence tagging can be exploited to mitigate the background due to $^{77}$Ge and its metastable state $^{77m}$Ge, which have been identified as dominant cosmogenic background in the search for neutrinoless double beta decay of $^{76}$Ge [282].

*   **Purification.** A number of high efficiency specific purification methods have been developed for different detectors in order to remove long-lived radio-isotopes, $^{39}$K, $^{85}$Kr, and $^{210}$Pb progeny. For organic liquid scintillators, distillation and water extraction have been shown to be very effective to reach radiopurity levels of the order of $10^{-5} \mu$Bq/kg or better in uranium and thorium, and $10^{-3} \mu$Bq/kg in $^{210}$Pb [283]. Distillation has been used to reduce $^{85}$Kr in xenon by a factor of $10^3$ [284]. Cryogenic distillation can be used to reduce the isotopic abundance of $^{39}$Ar in argon extracted from underground with a 350 m column in the ARIA project [285]. For semiconductors [286] and scintillating crystals [287] zone-refining followed by single-crystal growth has been exploited to remove impurities at the cost of a small fraction of material kept after the purification.

*   **Direct isotope tagging.** In neutrinoless double beta decay in order to explore half-lives greater than $10^{28}$ y one needs an almost background-free detector. The most robust signature is the identification of the daughter atom in the decay: for example, for $^{136}$Xe, the $^{136}$Ba$^{2+}$. Important step forwards to establish a valid and promising method for this tagging have been recently made [288,289].

### References

[1] A. Aab et al. (Pierre Auger), Earth and Space Science 7, 4, e2019EA000582 (2020), URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019EA000582.

[2] R. M. Baltrusaitis et al., Nucl. Instrum. Meth. A240, 410 (1985).

[3] D. J. Bird et al. (HiRes), Astrophys. J. 424, 491 (1994).

[4] T. Abu-Zayyad et al., Nucl. Instrum. Meth. A450, 253 (2000).

[5] A. Watson, European Phys. J. Web of Confs. 210, 00001 (2019), [arXiv:1901.06676].

[6] P. Sokolsky, European Phys. J. Web of Confs. 283, 01003 (2023).

[7] H. Tokuno et al., Nucl. Instrum. Meth. A676, 54 (2012), [arXiv:1201.0002].

[8] J. Abraham et al. (Pierre Auger), Nucl. Instrum. Meth. A620, 227 (2010), [arXiv:0907.4282].
[9] A. C. et al., Astroparticle Physics 149, 102819 (2023), ISSN 0927-6505, URL https://www.sciencedirect.com/science/article/pii/S0927650523000051.

[10] M. M. Malacari et al. (FAST), Astroparticle Physics 119, 102430 (2020), ISSN 0927-6505, URL https://www.sciencedirect.com/science/article/pii/S0927650520300037.

[11] F. Arqueros, J. R. Hoerandel and B. Keilhauer, Nucl. Instrum. Meth. A597, 23 (2008), [arXiv:0807.3844].

[12] F. Arqueros, J. R. Hoerandel and B. Keilhauer, Nucl. Instrum. Meth. A597, 1 (2008), [arXiv:0807.3760].

[13] J. Rosado, F. Blanco and F. Arqueros, Astropart. Phys. 34, 164 (2010), [arXiv:1004.3971].

[14] M. Ave et al. (AIRFLY), Astropart. Phys. 28, 41 (2007), [arXiv:astro-ph/0703132].

[15] J. Rosado, F. Blanco and F. Arqueros, Astropart. Phys. 55, 51 (2014), [arXiv:1401.4310].

[16] J. H. Boyer et al., Nucl. Instrum. Meth. A482, 457 (2002).

[17] J. T. Brack et al., Astropart. Phys. 20, 653 (2004).

[18] B. Fick et al., JINST 1, 11, P11003 (2006).

[19] J. Abraham et al. (Pierre Auger), Astropart. Phys. 33, 108 (2010), [arXiv:1002.0366].

[20] P. Abreu et al. (Pierre Auger), Astropart. Phys. 34, 368 (2011), [arXiv:1010.6162].

[21] R. Thalman et al. Journal of Quantitative Spectroscopy and Radiative Transfer, 147, 171 (2014), Erratum-ibid. 189, 281 (2017).

[22] J. Abraham et al. [Pierre Auger Collab.], Nucl. Instrum. Methods A789, 172 (2015).

[23] M. Unger et al., Nucl. Instrum. Meth. A588, 433 (2008), [arXiv:0801.4309].

[24] T.K. Gaisser and A.M. Hillas, Proc. 15th Int. Cosmic Ray Conf. Bulgarska Akademiia na Naukite, Conf. Papers 8, 353 (1978), (archived at http://adsabs.harvard.edu/abs/1977ICRC...8..353G).

[25] A. e. a. T. P. A. C. Aab (The Pierre Auger Collaboration), Phys. Rev. D 100, 082003 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.082003.

[26] A. Abdul Halim et al. (Pierre Auger), "Astrophys. J. Suppl." 264, 2, 50 (2023).

[27] A. H. et al. (Pierre Auger Collaboration), Phys. Rev. Lett. 132, 021001 (2024), URL https://link.aps.org/doi/10.1103/PhysRevLett.132.021001.

[28] P. A. Klimov et al., Space Sci. Rev. 212, 3-4, 1687 (2017), [arXiv:1706.04976].

[29] B. A. e. a. Khrenov, Cosmic Research 58, 5, 317 (2020).

[30] G. Abdellaoui et al. (JEM-EUSO), Journal of Instrumentation 13, 5 (2018), ISSN 17480221.

[31] G. Abdellaoui et al. (The JEM-EUSO Collaboration), Astroparticle Physics 154, 102891 (2024).

[32] J. Eser et al. (The JEM-EUSO Collaboration), Proceedings of Science 444, 397 (2023).

[33] F. Capel et al., Advances in Space Research 62, 2954 (2018).

[34] A. O. et al., Journal of Cosmology and Astroparticle Physics 2021, 06, 007 (2021), URL https://doi.org/10.1088/1475-7516/2021/06/007.

[35] J. Holder et al., AIP Conf. Proc. 1085, 657 (2009), [arXiv:0810.0474].

[36] F. Aharonian et al. (H.E.S.S.), Astron. Astrophys. 457, 899 (2006), [arXiv:astro-ph/0607333].

[37] J. Albert et al. (MAGIC), Astrophys. J. 674, 1037 (2008), [arXiv:0705.3244].

[38] T. C. Weekes et al., Astrophys. J. 342, 379 (1989).

[39] A. M. Hillas et al., Astrophys. J. 503, 744 (1998).

[40] http://tevcat.uchicago.edu/.

[41] F. Aharonian et al. (H.E.S.S.), Astrophys. J. 636, 777 (2006), [arXiv:astro-ph/0510397].

[42] M. de Naurois and D. Mazin, Comptes Rendus Physique 16, 610 (2015), [arXiv:1511.00463].

[43] A. Mitchell, in "37th International Cosmic Ray Conference," 46 (2022), [arXiv:2109.13753].

[44] A. M. Hillas, Astropart. Phys. 43, 19 (2013).

[45] M. Tluczykont et al., in "The Sixteenth Marcel Grossmann Meeting. On Recent Developments in Theoretical and Experimental General Relativity," 3324-3342 (2023).

[46] Z. Cao et al., ChA&A 43, 4, 457 (2019).

[47] B. S. Acharya et al. (CTA Consortium), Astropart. Phys. 43, 3 (2013).

[48] C.-L. Project et al., arXiv e-prints arXiv:2306.12960 (2023), [arXiv:2306.12960].

[49] Cherenkov Telescope Array Consortium et al., Science with the Cherenkov Telescope Array (2019).

[50] A. Bernstein et al., Report on the Depth Requirements for a Massive Detector at Homestake (2009), [arXiv:0907.4183].

[51] A. Suzuki, Eur. Phys. J. C 74, 10, 3094 (2014), [arXiv:1409.4515].

[52] G. Alimonti et al. (Borexino), Nucl. Instrum. Meth. A 600, 568 (2009), [arXiv:0806.2400].

[53] Y. Ashie et al. (Super-Kamiokande), Phys. Rev. D71, 112005 (2005), [hep-ex/0501064].

[54] S. Kasuga et al., Phys. Lett. B374, 238 (1996).

[55] M. Shiozawa (Super-Kamiokande), Nucl. Instrum. Meth. A433, 240 (1999).

[56] K. Abe et al. (Super-Kamiokande), Phys. Rev. D83, 052010 (2011), [arXiv:1010.0118].

[57] J. F. Beacom and M. R. Vagins, Phys. Rev. Lett. 93, 171101 (2004), [hep-ph/0309300].

[58] J. Boger et al. (SNO), Nucl. Instrum. Meth. A449, 172 (2000), [arXiv:nucl-ex/9910016].

[59] T. K. Gaisser, F. Halzen and T. Stanev, Phys. Rept. 258, 173 (1995), [Erratum: Phys. Rept.271,355(1996)], [hep-ph/9410384].

[60] J.G. Learned and K. Mannheim, Ann. Rev. Nucl. and Part. Sci. 50, 679 (2000).

[61] U. F. Katz and C. Spiering, Prog. Part. Nucl. Phys. 67, 651 (2012), [arXiv:1111.0507].

[62] M. G. Aartsen et al. (IceCube), J. Phys. G44, 5, 054006 (2017), [arXiv:1607.02671].

[63] S. Adrián-Martínez et al. (KM3NeT), J. Phys. G43, 8, 084001 (2016), [arXiv:1601.07459].

[64] A. Gazizov and M. P. Kowalski, Comput. Phys. Commun. 172, 203 (2005), [arXiv:astro-ph/0406439].

[65] A. D. Avrorin et al., Phys. Part. Nucl. 46, 2, 211 (2015).

[66] M. G. Aartsen et al. (IceCube), JINST 12, 03, P03012 (2017), [arXiv:1612.05093].

[67] M. G. Aartsen et al. (IceCube) (2014), [arXiv:1412.5106].

[68] M. G. Aartsen et al. (IceCube-Gen2), J. Phys. G 48, 6, 060501 (2021), [arXiv:2008.04323].

[69] M. Agostini et al. (P-ONE), Nature Astron. 4, 10, 913 (2020), [arXiv:2005.09493].

[70] Z. P. Ye et al. (2022), [arXiv:2207.04519].

[71] F. Hu, PoS ICRC2023, 1203 (2023).

[72] T. Huege, Phys. Rept. 620, 1 (2016), [arXiv:1601.07426].

[73] F. G. Schroder, Prog. Part. Nucl. Phys. 93, 1 (2017), [arXiv:1607.08781].

[74] G. A. Askar'yan, Sov. Phys. JETP 14, 2, 441 (1962), [Zh. Eksp. Teor. Fiz.41,616(1961)].
[75] G.A. Askaryan, Sov. Phys. JETP **21**, 658 (1965).

[76] E. Zas, F. Halzen and T. Stanev, Phys. Rev. D **45**, 362 (1992).

[77] C. W. James et al., Phys. Rev. E**84**, 056602 (2011), [arXiv:1007.4146].

[78] C. W. James, Phys. Rev. D **105**, 2, 023014 (2022), [arXiv:2201.01298].

[79] J. Alvarez-Muniz, R. A. Vazquez and E. Zas, Phys. Rev. D**62**, 063001 (2000), [arXiv:astro-ph/0003315].

[80] C. Glaser et al., Eur. Phys. J. C **80**, 77 (2020), [arXiv:1906.01670].

[81] D. Saltzberg et al., Phys. Rev. Lett. **86**, 2802 (2001), [hep-ex/0011001].

[82] O. Scholten et al., J. Phys. Conf. Ser. **81**, 012004 (2007).

[83] J. Alvarez-Muñiz et al., Phys. Rev. D **74**, 023007 (2006), URL https://link.aps.org/doi/10.1103/PhysRevD.74.023007.

[84] L. Gerhardt and S. R. Klein, Phys. Rev. D**82**, 074017 (2010), [arXiv:1007.0039].

[85] J. Alvarez-Muniz, R. A. Vazquez and E. Zas, Phys. Rev. D**61**, 023001 (2000), [arXiv:astro-ph/9901278].

[86] D. García-Fernández, A. Nelles and C. Glaser, Phys. Rev. D **102**, 8, 083011 (2020), [arXiv:2003.13442].

[87] K. Greisen, Phys. Rev. Lett. **16**, 748 (1966).

[88] G. T. Zatsepin and V. A. Kuzmin, JETP Lett. **4**, 78 (1966), [Pisma Zh. Eksp. Teor. Fiz.4,114(1966)].

[89] V. de Souza (Pierre Auger, Telescope Array), PoS ICRC2017, 522 (2018).

[90] R. Abbasi et al. (IceCube), Astropart. Phys. **34**, 382 (2011), [arXiv:1004.1694].

[91] P. Allison et al. (ARA), Astropart. Phys. **108**, 63 (2019), [arXiv:1712.03301].

[92] A. Anker et al. (ARIANNA), JCAP **1911**, 030 (2019), [arXiv:1909.02677].

[93] J. A. Aguilar et al., Eur. Phys. J. C **82**, 2, 147 (2022), [arXiv:2107.02604].

[94] I. Plaisier, S. Bouma and A. Nelles, Eur. Phys. J. C **83**, 5, 443 (2023), [arXiv:2302.00054].

[95] A. G. Vieregg, K. Bechtol and A. Romero-Wolf, JCAP **1602**, 02, 005 (2016), [arXiv:1504.08006].

[96] J. A. Aguilar et al., Astropart. Phys. **145**, 102790 (2023), [arXiv:2103.06079].

[97] S. W. Barwick et al. (ARIANNA), Astropart. Phys. **90**, 50 (2017), [arXiv:1612.04473].

[98] M. G. Aartsen et al. (IceCube), Phys. Rev. D**98**, 6, 062003 (2018), [arXiv:1807.01820].

[99] A. Aab et al. (Pierre Auger), JCAP **1910**, 10, 022 (2019), [arXiv:1906.07422].

[100] J. D. Bray et al., Phys. Rev. D**91**, 6, 063002 (2015), [arXiv:1502.03313].

[101] O. Scholten et al., Phys. Rev. Lett. **103**, 191301 (2009), [arXiv:0910.4745].

[102] P. W. Gorham et al. (ANITA), Phys. Rev. D**99**, 12, 122001 (2019), [arXiv:1902.04005].

[103] P. Allison et al. (ARA), Phys. Rev. D **102**, 4, 043021 (2020), [arXiv:1912.00987].

[104] A. Anker et al. (ARIANNA), JCAP **03**, 053 (2020), [arXiv:1909.00840].

[105] T. Winchen et al., J. Phys. Conf. Ser. **1181**, 1, 012077 (2019), [arXiv:1903.08472].

[106] R. Abbasi et al. (IceCube), Astrophys. J. **928**, 1, 50 (2022), [arXiv:2111.10299].

[107] A. van Vliet, R. Alves Batista and J. R. Hörandel, Phys. Rev. D **100**, 2, 021302 (2019), [arXiv:1901.01899].

[108] S. Barwick et al., Journal of Glaciology **51**, 173, 231–238 (2005).

[109] J.A. Dowdeswell and S. Evans, Rept. on Prog. in Phys. **67**, 1821 (2004).

[110] S. W. Barwick et al., JCAP **1807**, 07, 055 (2018), [arXiv:1804.10430].

[111] C. Deaconu et al., Phys. Rev. D**98**, 4, 043010 (2018), [arXiv:1805.12576].

[112] P. W. Gorham et al. (ANITA), Phys. Rev. Lett. **103**, 051103 (2009), [arXiv:0812.2715].

[113] P. W. Gorham et al. (ANITA), Phys. Rev. Lett. **121**, 16, 161102 (2018), [arXiv:1803.05088].

[114] P. W. Gorham et al. (ANITA), Phys. Rev. Lett. **126**, 7, 071103 (2021), [arXiv:2008.05690].

[115] I. Kravchenko et al. (RICE), Phys. Rev. D**73**, 082002 (2006), [arXiv:astro-ph/0601148].

[116] P. Allison et al. (ARA), Phys. Rev. D**93**, 8, 082003 (2016), [arXiv:1507.08991].

[117] J. Avva et al., Nucl. Instrum. Meth. **A869**, 46 (2017), [arXiv:1605.03525].

[118] J. A. Aguilar et al. (RNO-G), JINST **16**, 03, P03025 (2021), [arXiv:2010.12279].

[119] R.D. Dagkesamanskii and I.M. Zheleznykh, Sov. Phys. JETP Lett. **50**, 233 (1989).

[120] G. R. Olhoeft and D. W. Strangway, Earth and Planetary Science Letters **24**, 3, 394 (1975).

[121] J. D. Bray, Astropart. Phys. **77**, 1 (2016), [arXiv:1601.02980].

[122] C. W. James et al., EPJ Web Conf. **135**, 04001 (2017), [arXiv:1704.05336].

[123] A. Aab et al. (Pierre Auger), Phys. Rev. Lett. **116**, 24, 241101 (2016), [arXiv:1605.02564].

[124] P. A. Bezyazeekov et al. (Tunka-Rex), Phys. Rev. D**97**, 12, 122004 (2018), [arXiv:1803.06862].

[125] S. Buitink et al. (LOFAR), Phys. Rev. D**90**, 8, 082003 (2014), [arXiv:1408.7001].

[126] K. Mulrey et al. (LOFAR), JCAP **11**, 017 (2020), [arXiv:2005.13441].

[127] D. Fargion, Astrophys. J. **570**, 909 (2002), [arXiv:astro-ph/0002453].

[128] X. Bertou et al., Astropart. Phys. **17**, 183 (2002), [arXiv:astro-ph/0104452].

[129] J. L. Feng et al., Phys. Rev. Lett. **88**, 161102 (2002), [hep-ph/0105067].

[130] Q. Abarr et al. (PUEO), JINST **16**, 08, P08035 (2021), [arXiv:2010.02892].

[131] J. Álvarez Muñiz et al. (GRAND), Sci. China Phys. Mech. Astron. **63**, 1, 219501 (2020), [arXiv:1810.09994].

[132] P. Abreu et al. (Pierre Auger), PoS ICRC2021, 262 (2021).

[133] D. Nygren, eConf C**740805**, 58 (1974).

[134] V. Chepel and H. Araujo, JINST **8**, R04001 (2013), [arXiv:1207.2292].

[135] E. Aprile and T. Doke, Rev. Mod. Phys. **82**, 2053 (2010), [arXiv:0910.4956].

[136] D. Gonzalez-Diaz, F. Monrabal and S. Murphy, Nucl. Instrum. Meth. A **878**, 200 (2018), [arXiv:1710.01018].

[137] E. Aprile et al., Noble Gas Detectors, Wiley (2008), ISBN 978-3-527-40597-8, 978-3-527-61002-0.

[138] A. I. Bolozdynya, Emission detectors, World Scientific (2010).

[139] D. S. Akerib et al. (LUX), Phys. Rev. D **97**, 11, 112002 (2018), [arXiv:1802.06162].

[140] B. Lenardo et al., IEEE Trans. Nucl. Sci. **62**, 6, 3387 (2015), [arXiv:1412.4417].
[141] W. F. Schmidt, IEEE Trans. Electric. Insul. **19**, 389 (1984).

[142] L. Baudis et al., Eur. Phys. J. C **83**, 8, 717 (2023), [arXiv:2303.13963].

[143] J. Albert et al. (EXO-200), Phys. Rev. C **89**, 1, 015502 (2014), [arXiv:1306.6106].

[144] A. Gando et al. (KamLAND-Zen), Phys. Rev. Lett. **122**, 19, 192501 (2019), [arXiv:1901.03871].

[145] E. Aprile et al. (XENON), Phys. Rev. C **106**, 2, 024328 (2022), [arXiv:2205.04158].

[146] J. B. Albert et al. (EXO-200), Phys. Rev. D **96**, 9, 092001 (2017), [arXiv:1704.05042].

[147] K. Abe et al. (XMASS), PTEP **2018**, 5, 053D03 (2018), [arXiv:1801.03251].

[148] T. Doke et al., Jap. J. Appl. Phys. **41**, 1538 (2002).

[149] C. E. Dahl, *The Physics of Background Discrimination in Liquid Xenon, and First Results from XENON10 in the Hunt for WIMP Dark Matter*, Ph.D. thesis, Princeton University (2009).

[150] L. Baudis, P. Sanchez-Lucas and K. Thieme, Eur. Phys. J. C **81**, 12, 1060 (2021), [arXiv:2109.07151].

[151] G. Anton et al. (EXO-200), Phys. Rev. C **101**, 6, 065501 (2020), [arXiv:1908.04128].

[152] J. J. Gomez-Cadenas, F. Monrabal Capilla and P. Ferrario, Front. in Phys. **7**, 51 (2019), [arXiv:1903.02435].

[153] M. S. J. Lindhard, V. Nielsen and P. Thomsen, Kgl. Danske Videnskab., Selskab. Mat. Fys. Medd **33**, 1 (1963).

[154] A. Hitachi and A. Mozumder (2019), [arXiv:1903.05815].

[155] P. Sorensen and C. E. Dahl, Phys. Rev. D **83**, 063501 (2011), [arXiv:1101.6080].

[156] M. Szydagis et al., JINST **6**, P10002 (2011), [arXiv:1106.1613].

[157] G. Pereira, C. Silva and V. N. Solovov (LZ), JINST **18**, 04, C04007 (2023).

[158] E. Aprile et al. (XENON), Eur. Phys. J. C **80**, 8, 785 (2020), [arXiv:2003.03825].

[159] G. Anton et al. (EXO-200), Phys. Rev. Lett. **123**, 16, 161802 (2019), [arXiv:1906.02723].

[160] D. Nygren, Nucl. Instrum. Meth. A **603**, 337 (2009).

[161] J. Renner et al. (NEXT), JHEP **10**, 230 (2019), [arXiv:1905.13110].

[162] P. Agnes et al. (DarkSide), Nucl. Instrum. Meth. A **904**, 23 (2018), [arXiv:1802.01427].

[163] P. Abratenko et al. (MicroBooNE), JINST **16**, 09, P09025 (2021), [arXiv:2104.06551].

[164] Y. Li et al., Nucl. Instrum. Meth. A **816**, 160 (2016), [arXiv:1508.07059].

[165] J. B. Albert et al. (EXO-200), Phys. Rev. C **95**, 2, 025502 (2017), [arXiv:1609.04467].

[166] M. Babicz et al., JINST **15**, 09, P09009 (2020), [arXiv:2002.09346].

[167] E. Aprile et al. (XENON), PTEP **2022**, 5, 053H01 (2022), [arXiv:2112.12231].

[168] P. Agnes et al. (DarkSide), Phys. Rev. D **98**, 10, 102006 (2018), [arXiv:1802.07198].

[169] E. Aprile et al. (XENON), Phys. Rev. Lett. **129**, 16, 161805 (2022), [arXiv:2207.11330].

[170] P. Agnes (DarkSide), EPJ Web Conf. **280**, 06003 (2023).

[171] D. S. Akerib et al. (LZ), Nucl. Instrum. Meth. A **953**, 163047 (2020), [arXiv:1910.09124].

[172] Y. Meng et al. (PandaX-4T), Phys. Rev. Lett. **127**, 26, 261802 (2021), [arXiv:2107.13438].

[173] E. Aprile et al. (XENON), JCAP **11**, 031 (2020), [arXiv:2007.08796].

[174] S. Kravitz et al., Eur. Phys. J. C **80**, 3, 262 (2020), [arXiv:1909.08730].

[175] J. Aalbers et al. (DARWIN), JCAP **1611**, 11, 017 (2016), [arXiv:1606.07001].

[176] J. Aalbers et al., J. Phys. G **50**, 1, 013001 (2023), [arXiv:2203.02309].

[177] X. Wang et al., JINST **18**, 05, P05028 (2023), [arXiv:2301.06044].

[178] P. Agnes et al. (DarkSide), Phys. Rev. Lett. **130**, 10, 101002 (2023), [arXiv:2207.11968].

[179] S. Li et al. (PandaX), Phys. Rev. Lett. **130**, 26, 261001 (2023), [arXiv:2212.10067].

[180] E. Aprile et al. (XENON), Phys. Rev. Lett. **129**, 16, 161805 (2022), [arXiv:2207.11330].

[181] J. Aalbers et al. (LZ) (2023), [arXiv:2307.15753].

[182] T. N. Thorpe (DarkSide-20k), SciPost Phys. Proc. **12**, 069 (2023), [arXiv:2210.00322].

[183] G. Plante et al., Eur. Phys. J. C **82**, 10, 860 (2022), [arXiv:2205.07336].

[184] A. Abdukerim et al. (PandaX) (2024), [arXiv:2402.03596].

[185] K. Stifter (LZ), J. Phys. Conf. Ser. **1468**, 1, 012016 (2020).

[186] R. Anthony-Petersen et al. (2023), [arXiv:2307.11877].

[187] B. von Krosigk et al., SciPost Phys. Proc. **12**, 016 (2023), [arXiv:2209.10950].

[188] R. E. Lanou, H. J. Maris and G. M. Seidel, Phys. Rev. Lett. **58**, 2498 (1987).

[189] F. D. Amaro et al., Nucl. Instrum. Meth. A **1054**, 168325 (2023), [arXiv:2306.04568].

[190] G. Anton et al. (EXO-200), Phys. Rev. Lett. **123**, 16, 161802 (2019), [arXiv:1906.02723].

[191] G. Adhikari et al. (nEXO), J. Phys. G **49**, 1, 015104 (2022), [arXiv:2106.16243].

[192] A. Simón et al. (NEXT), JHEP **07**, 146 (2021), [arXiv:2102.11931].

[193] P. Ferrario et al. (NEXT), JHEP **10**, 052 (2019), [arXiv:1905.13141].

[194] P. Novella et al. (NEXT), Phys. Rev. C **105**, 5, 055501 (2022), [arXiv:2111.11091].

[195] C. Adams et al. (NEXT), JHEP **2021**, 08, 164 (2021), [arXiv:2005.06467].

[196] A. D. McDonald et al., Phys. Rev. Lett. **120**, 13, 132504 (2018), [arXiv:1711.04782].

[197] P. Herrero-Gómez et al. (NEXT), Nature Commun. **13**, 1, 7741 (2022), [arXiv:2201.09099].

[198] H. Rasiwala et al. (nEXO Ba-tagging group), Nucl. Instrum. Meth. B **541**, 298 (2023), [arXiv:2303.04698].

[199] T. Li et al., JHEP **05**, 200 (2023), [arXiv:2211.14992].

[200] E. Aprile et al. (XENON), Nature **568**, 7753, 532 (2019), [arXiv:1904.11002].

[201] G. Martínez-Lema et al. (NEXT), JHEP **21**, 203 (2020), [arXiv:2006.07320].

[202] C. Wittweg et al., Eur. Phys. J. C **80**, 12, 1161 (2020), [arXiv:2002.04239].

[203] L. Si et al. (PandaX), Research **2022**, 9798721 (2022), [arXiv:2205.12809].

[204] D. S. Akerib et al. (LZ), Phys. Rev. C **102**, 1, 014602 (2020), [arXiv:1912.04248].

[205] D. S. Akerib et al. (LZ), Phys. Rev. C **104**, 6, 065501 (2021), [arXiv:2104.13374].

[206] F. Agostini et al. (DARWIN), Eur. Phys. J. C **80**, 9, 808 (2020), [arXiv:2003.13407].

[207] C. Enss, editor, Cryogenic particle detection, volume 99 of Topics in applied physics, Springer, Berlin, Germany (2005).
[208] K. Pretzl, *Cryogenic Detectors* (2020).

[209] Proc. of the Low Temperature Detectors for Neutrinos and Dark Matter, Low Temperature Detectors for Neutrinos and Dark Matter (start 1987).

[210] G. H. Rieke, *Detection of Light*, Cambridge University Press, 3rd edition (2021).

[211] H. Kraus, Superconductor Science and Technology **9**, 10, 827 (1996), URL https://doi.org/10.1088/0953-2048/9/10/001.

[212] A. V. Dixit et al., Phys. Rev. Lett. **126**, 14, 141302 (2021), [arXiv:2008.12231].

[213] J. Zmuidzinas, Annual Review of Condensed Matter Physics **3**, 1, 169 (2012), URL https://doi.org/10.1146/annurev-conmathphys-020911-125022.

[214] S. Doyle et al., Journal of Low Temperature Physics **151**, 1, 530 (2008), URL https://doi.org/10.1007/s10909-007-9685-2.

[215] K. Irwin and G. Hilton, Transition-Edge Sensors, 63–150, Springer Berlin Heidelberg, Berlin, Heidelberg (2005), ISBN 978-3-540-31478-3, URL https://doi.org/10.1007/10933596_3.

[216] F. Paolucci and F. Giazotto, Instruments **5**, 2 (2021), ISSN 2410-390X, URL https://www.mdpi.com/2410-390X/5/2/14.

[217] S. B. Bandler et al., J. Low Temp. Phys. **93**, 709 (1993).

[218] D. Prele, Journal of Instrumentation **10**, 08, C08015 (2015), URL https://doi.org/10.1088/1748-0221/10/08/c08015.

[219] E. Fiorini and T. Niinikoski, Nucl. Instrum. Meth. A **224**, 83 (1984).

[220] S. H. Moseley, J. C. Mather and D. McCammon, Thermal detectors as x-ray spectrometers (1984).

[221] V. Barger and D. Cline (1985), telemark, 1984.

[222] A. Blasi et al. (1985), I.N.F.N./BE-85/2, internal report.

[223] D. Q. Adams et al. (CUORE) (2021), [arXiv:2104.06906].

[224] W. R. Armstrong et al. (CUPID) (2019), [arXiv:1907.09376].

[225] O. Azzolini et al. (CUPID), Phys. Rev. Lett. **123**, 3, 032501 (2019), [arXiv:1906.05001].

[226] E. Armengaud et al. (CUPID), Phys. Rev. Lett. **126**, 18, 181802 (2021), [arXiv:2011.13243].

[227] V. Alenkov et al., Eur. Phys. J. C **79**, 9, 791 (2019), [arXiv:1903.09483].

[228] D. McCammon et al. (1984).

[229] P. Filianin et al., Phys. Rev. Lett. **127**, 7, 072502 (2021), [arXiv:2108.07039].

[230] D. Pergolesi et al., Nucl. Instrum. Meth. A **559**, 349 (2006).

[231] M. Sisti et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **520**, 1, 125 (2004), ISSN 0168-9002, proceedings of the 10th International Workshop on Low Temperature Detectors, URL https://www.sciencedirect.com/science/article/pii/S0168900203031814.

[232] E. Ferri et al., Phys. Procedia **61**, 227 (2015).

[233] C. Velte, Measurement of a high energy resolution and high-statistics 163Ho electron capture spectrum for the ECHO experiment., Ph.D. thesis, U. Heidelberg (main) (2020).

[234] L. Gastaldo et al., Eur. Phys. J. ST **226**, 8, 1623 (2017).

[235] B. Alpert et al., Eur. Phys. J. C **75**, 3, 112 (2015), [arXiv:1412.5060].

[236] M. P. Croce et al., J. Low Temp. Phys. **184**, 3-4, 958 (2016), [arXiv:1510.03874].

[237] R. Ardito et al. (2005), [hep-ex/0501010].

[238] P. Luke et al., Nucl. Instrum. Meth. A **289**, 406 (1990).

[239] B. Neganov et al., J Low Temp Phys **93**, 417-422 (1993), URL https://doi.org/10.1007/BF00693454.

[240] I. Alkhatib et al. (SuperCDMS), Phys. Rev. Lett. **127**, 8, 081802 (2021), [arXiv:2011.09183].

[241] Q. Arnaud et al. (EDELWEISS), Phys. Rev. Lett. **125**, 14, 141301 (2020), [arXiv:2003.01046].

[242] Snowmass2021-Letter of Interest The TESSERACT Dark Matter Project (2020).

[243] T. Braine et al. (ADMX), Phys. Rev. Lett. **124**, 10, 101303 (2020), [arXiv:1910.08638].

[244] Snowmass2021-Letter of Interest Opening the Terahertz Axion Window (2020).

[245] C. Bellenghi et al., Eur. Phys. J. C **79**, 9, 727 (2019), [arXiv:1905.10611].

[246] I. Colantoni et al., J. Low Temp. Phys. **199**, 3-4, 593 (2020).

[247] G. Heusser, Ann. Rev. Nucl. Part. Sci. **45**, 543 (1995).

[248] D. S. Akerib et al., The European Physical Journal C **80**, 11 (2020), URL https://doi.org/10.1140/2Fepjc/2Fs10052-020-8420-x.

[249] D. S. Leonard et al., Nucl. Instrum. Meth. A **871**, 169 (2017), [arXiv:1703.10799].

[250] N. Abgrall et al. (Majorana), Nucl. Instrum. Meth. A **823**, 83 (2016), [arXiv:1603.08483].

[251] E. Aprile et al. (XENON), Eur. Phys. J. C **82**, 7, 599 (2022), [arXiv:2112.05629].

[252] M. Laubenstein, International Journal of Modern Physics A **32**, 30, 1743002 (2017), URL https://doi.org/10.1142/s0217751x17430023.

[253] P. Scovell et al., Astroparticle Physics **97**, 160 (2018), URL https://doi.org/10.1016/j.astropartphys.2017.11.006.

[254] G. Heusser et al., The European Physical Journal C **75**, 11 (2015), URL https://doi.org/10.1140/2Fepjc/2Fs10052-015-3704-2.

[255] N. Jakubowski, Analytical and Bioanalytical Chemistry **392**, 5, 775 (2008), URL https://doi.org/10.1007/2Fs00216-008-2374-4.

[256] M. Clemenza, J. Radioanal. Nucl. Chem. **318**, 3, 1765 (2018).

[257] G. Heusser et al., Applied Radiation and Isotopes **52**, 3, 691 (2000), URL https://doi.org/10.1016/2Fs0969-8043/2899/2900231-6.

[258] M. Wojcik, G. Zuzel and H. Simgen, International Journal of Modern Physics A **32**, 30, 1743004 (2017), URL https://doi.org/10.1142/2Fs0217751x17430047.

[259] Y. Takeuchi et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **421**, 1-2, 334 (1999), URL https://doi.org/10.1016/2Fs0168-9002/2898/2901204-2.

[260] J. Kiko, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **460**, 2-3, 272 (2001), URL https://doi.org/10.1016/2Fs0168-9002/2800/2901082-2.

[261] E. Hoppe et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **579**, 1, 486 (2007), URL https://doi.org/10.1016/2Fj.nima.2007.04.101.

[262] R. Bunker et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **967**, 163870 (2020), URL https://doi.org/10.1016/j.nima.2020.163870.

[263] S. Bruenner et al., The European Physical Journal C **81**, 4 (2021), URL https://doi.org/10.1140/epjc/s10052-021-09047-2.
[264] M. L. di Vacri et al., Nucl. Instrum. Meth. A **994**, 165051 (2021), [arXiv:2006.12746].

[265] H. Loosli, Earth and Planetary Science Letters **63**, 1, 51 (1983), URL https://doi.org/10.1016%2F0012-821x%2883%2990021-3.

[266] D. Acosta-Kane et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **587**, 1, 46 (2008), URL https://doi.org/10.1016%2Fj.nima.2007.12.032.

[267] P. Agnes et al. (DarkSide), Phys. Rev. D **93**, 8, 081101 (2016), [Addendum: Phys.Rev.D 95, 069901 (2017)], [arXiv:1510.00702].

[268] H. O. Back et al., in "Snowmass 2021," (2022), [arXiv:2203.09734].

[269] I. Bandac et al., Appl. Radiat. Isot. **126**, 127 (2017).

[270] G. Ranucci, in "IEEE Symposium Conference Record Nuclear Science 2004," volume 2, 804-809 Vol. 2 (2004).

[271] P. Adhikari et al., The European Physical Journal C **80**, 4 (2020), URL https://doi.org/10.1140%2Fepjc%2Fs10052-020-7789-x.

[272] D. Budjáš et al., Journal of Instrumentation **4**, 10, P10007 (2009), URL https://doi.org/10.1088/1748-0221/4/10/p10007.

[273] M. Agostini et al. (GERDA), Phys. Rev. Lett. **125**, 25, 252502 (2020), [arXiv:2009.06079].

[274] A. Simón et al. (NEXT), JHEP **21**, 146 (2020), [arXiv:2102.11931].

[275] R. Agnese et al. (SuperCDMS Collaboration), Phys. Rev. Lett. **120**, 061802 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.061802.

[276] W. Westphal et al., Nucl. Instrum. Meth. A **559**, 372 (2006).

[277] S. Westdale, E. Shields and F. Calaprice, Astroparticle Physics **79**, 10 (2016), URL https://doi.org/10.1016%2Fj.astropartphys.2016.01.005.

[278] Marco Selvi, The XENONnT Neutron Veto (2023), URL https://indico.cern.ch/event/1188759/.

[279] D. Akimov et al. (COHERENT), Science **357**, 6356, 1123 (2017), [arXiv:1708.01294].

[280] D. Akimov et al., in "Snowmass 2021," (2022), [arXiv:2204.04575].

[281] H. Back et al. (Borexino Collaboration), Phys. Rev. C **74**, 045805 (2006), URL https://link.aps.org/doi/10.1103/PhysRevC.74.045805.

[282] C. Wiesinger, L. Pandola and S. Schönert, Eur. Phys. J. C **78**, 7, 597 (2018), [arXiv:1802.05040].

[283] J. Benziger et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment **587**, 2-3, 277 (2008), URL https://doi.org/10.1016/j.nima.2007.12.043.

[284] K. Abe et al., Nucl. Instrum. Meth. A **716**, 78 (2013), [arXiv:1301.2815].

[285] E. Aaron et al. (DarkSide-20k), Eur. Phys. J. C **83**, 5, 453 (2023), [arXiv:2301.09639].

[286] K.-P. Gradwohl et al., Journal of Instrumentation **15**, 12, P12010 (2020), URL https://doi.org/10.1088/1748-0221/15/12/p12010.

[287] B. Suerfu, F. Calaprice and M. Souza, Physical Review Applied **16**, 1 (2021), URL https://doi.org/10.1103/physrevapplied.16.014060.

[288] I. Rivilla et al., Nature **583**, 7814, 48 (2020), URL https://doi.org/10.1038/s41586-020-2431-5.

[289] A. McDonald et al., Physical Review Letters **120**, 13 (2018), URL https://doi.org/10.1103/physrevlett.120.132504.
# 37. Radioactivity and Radiation Protection

Revised July 2023 by S. Roesler (CERN) and M. Silari.

## 37.1 Definitions [1–3]

It would be desirable if legal protection limits could be expressed in directly measurable physical quantities. However, this does not allow quantifying biological effects of the exposure of the human body and its detriment to ionizing radiation.

For this reason, dose limits are expressed in terms of so-called protection quantities which, although calculable, are not measurable. Protection quantities are used to quantify the extent of exposure of the human body to ionizing radiation from both whole and partial body external irradiation and from intakes of radionuclides.

In order to demonstrate compliance with dose limits, so-called operational quantities are typically used, which aim at providing conservative estimates of protection quantities. Often radiation protection detectors used for individual and area monitoring are calibrated in terms of operational quantities and, thus, these quantities become “measurable”.

### 37.1.1 Physical quantities

* **Fluence, $\Phi$** (unit: $1/m^2$): The fluence is the quotient of the sum of the particle track lengths $dl$ in the volume $dV$

$$ \Phi = dl/dV . $$

It can also be expressed in terms of number of particles $dN$ incident upon a small sphere of cross-sectional area $da$

$$ \Phi = dN/da . $$

* **Absorbed dose, $D$** (unit: gray, 1 Gy=1 J/kg=100 rad): The absorbed dose is the energy imparted by ionizing radiation in a volume element of a specified material divided by the mass of this volume element.

* **Kerma, $K$** (unit: gray): Kerma is the sum of the initial kinetic energies of all charged particles set in motion by indirectly ionizing radiation in a volume element of the specified material divided by the mass of this volume element.

* **Linear energy transfer, $L$ or LET** (unit: J/m, often given in keV/µm, 1 keV/µm ≈ 1.602 × 10⁻¹⁰ J/m): The linear energy transfer is the mean energy, $dE$, lost by a charged particle owing to collisions with electrons in traversing a distance $dl$ in matter. Low-LET radiation: X rays and gamma rays (accompanied by charged particles due to interactions with the surrounding medium) or light charged particles such as electrons that produce sparse ionizing events far apart at a molecular scale ($L < 10$ keV/µm). High-LET radiation: neutrons and heavy charged particles that produce ionizing events densely spaced at a molecular scale ($L > 10$ keV/µm). While the above LET definition refers to electronic stopping power only, at low energy nuclear stopping power could be a significant fraction of the total stopping power.

* **Activity, $A$** (unit: Becquerel, 1 Bq=1/s=27 pCi): Activity is the expectation value of the number of nuclear decays occurring in a given quantity of material per unit time.

### 37.1.2 Protection quantities

* **Organ absorbed dose, $D_T$** (unit: gray): The mean absorbed dose in an organ or tissue $T$ of mass $m_T$ is defined as

$$ D_T = \frac{1}{m_T} \int_{m_T} Ddm . $$

* **Equivalent dose, $H_T$** (unit: sievert, 1 Sv=1 J/kg=100 rem):

The equivalent dose $H_T$ in an organ or tissue $T$ is equal to the sum of the absorbed doses $D_{T,R}$ in the organ or tissue caused by different radiation types $R$ weighted with so-called radiation weighting factors $w_R$:

$$ H_T = \sum_R w_R \times D_{T,R} . $$

It expresses long-term risks (primarily cancer and leukemia) from low-level chronic exposure. The values for $w_R$ recommended by ICRP [2] are given in Table 37.1.

<table><caption>Table 37.1: Radiation weighting factors, $w_R$.</caption><thead><tr><th>Radiation type</th><th>$w_R$</th></tr></thead><tbody><tr><td>Photons, electrons and muons</td><td>1</td></tr><tr><td>Neutrons, $E_n < 1$ MeV</td><td>2.5 + 18.2 × exp[-(ln $E_n$)²/6]</td></tr><tr><td>1 MeV ≤ $E_n$ ≤ 50 MeV</td><td>5.0 + 17.0 × exp[-(ln(2$E_n$))²/6]</td></tr><tr><td>$E_n$ > 50 MeV</td><td>2.5 + 3.25 × exp[-(ln(0.04$E_n$))²/6]</td></tr><tr><td>Protons and charged pions</td><td>2</td></tr><tr><td>Alpha particles, fission</td><td></td></tr><tr><td>fragments, heavy ions</td><td>20</td></tr></tbody></table>

* **Effective dose, $E$** (unit: sievert): The sum of the equivalent doses, weighted by the tissue weighting factors $w_T$ ($\sum_T w_T = 1$) of several organs and tissues $T$ of the body that are considered to be most sensitive [2], is called “effective dose”:

$$ E = \sum_T w_T \times H_T . \qquad (37.1) $$

### 37.1.3 Operational quantities

* **Dose equivalent, $H$** (unit: sievert): The dose equivalent at a point in tissue is given by:

$$ H = D \times Q \qquad (37.2) $$

where $D$ is the absorbed dose and $Q$ is the quality factor at that point. The quality factor at a point in tissue, is given by:

$$ Q = \frac{1}{D} \int_{L=0}^{\infty} Q(L) D_L dL $$

where $D_L$ is the distribution of $D$ in unrestricted linear energy transfer $L$ at the point of interest, and $Q(L)$ is the quality factor as a function of $L$. The integration is to be performed over $D_L$, due to all charged particles, excluding their secondary electrons.

* **Ambient dose equivalent, $H^*(10)$** (unit: sievert): The dose equivalent at a point in a radiation field that would be produced by the corresponding expanded and aligned field in a 30 cm diameter sphere of unit density tissue (so-called ICRU sphere/tissue with a mass composition of 76.2% oxygen, 11.1% carbon, 10.1% hydrogen and 2.6% nitrogen) at a depth of 10 mm on the radius vector opposing the direction of the aligned field. Ambient dose equivalent is the operational quantity for area monitoring.

* **Personal dose equivalent, $H_p(d)$** (unit: sievert): The dose equivalent in ICRU tissue at an appropriate depth, $d$, below a specified point on the human body. The specified point is normally taken to be where the individual dosimeter is worn. For the assessment of effective dose, $H_p(10)$ with a depth $d = 10$ mm is chosen, and for the assessment of the dose to the skin and to the hands and feet the personal dose equivalent, $H_p(0.07)$, with a depth $d = 0.07$ mm, is used. Personal dose equivalent is the operational quantity for individual monitoring.

### 37.1.4 Dose conversion coefficients

Dose conversion coefficients allow direct calculation of protection or operational quantities from particle fluence and are functions of particle type, energy and irradiation configuration. The most common coefficients are those for effective dose and ambient dose equivalent. The former are based on simulations in which the dose to organs of anthropomorphic phantoms is calculated for approximate actual conditions of exposure, such as irradiation of the front of the body (antero-posterior irradiation) or isotropic irradiation.

Conversion coefficients from fluence to effective dose are given for anterior-posterior irradiation and various particle types in Fig. 37.1 [4]. For example, the effective dose from an anterior-posterior irradiation in a field of 1-MeV neutrons with a fluence of 1 neutron per cm² is about 290 pSv. In Monte Carlo simulations such coefficients allow multiplication with fluence at scoring time such that effective dose to a human body at the considered location is directly obtained.
<figure><img src="image_251.png" /><figcaption>**Figure 37.1:** Fluence to effective dose conversion coefficients for anterior-posterior irradiation and various particle types [4].</figcaption></figure>

## 37.2 Radiation levels [5]

* **Natural background radiation:**

On a worldwide average, the annual whole-body dose equivalent due to all sources of natural background radiation ranges from 1.0 to 13 mSv (0.1–1.3 rem) with an annual average of 2.4 mSv [6]. In certain areas values up to 50 mSv (5 rem) have been measured. A large fraction (typically more than 50%) originates from inhaled natural radioactivity, mostly radon and radon daughters [7]. The latter can vary by more than one order of magnitude: it is 0.1–0.2 mSv in open areas, 2 mSv on average in a house and more than 20 mSv in poorly ventilated mines.

* **Cosmic ray background radiation:** At sea level, the whole-body dose equivalent due to cosmic ray background radiation is dominated by muons; at higher altitudes also nucleons contribute. Dose equivalent rates range from less than 0.1 µSv/h at sea level to a few µSv/h at aircraft altitudes. Details on cosmic ray fluence levels are given in the Cosmic Rays section (Sec. 30 of this Review).

## 37.3 Health effects of ionizing radiation

Radiation can cause two types of health effects, deterministic and stochastic:

* **Deterministic effects** are tissue reactions which cause injury to a population of cells if a given threshold of absorbed dose is exceeded. The severity of the reaction increases with dose. The quantity in use for tissue reactions is the absorbed dose, *D*. When particles other than photons and electrons (low-LET radiation) are involved, a Relative Biological Effectiveness (RBE)-weighted dose may be used. The RBE of a given radiation is the reciprocal of the ratio of the absorbed dose of that radiation to the absorbed dose of a reference radiation (usually X rays) required to produce the same degree of biological effect. It is a complex quantity that depends on many factors such as cell type, dose rate, fractionation, etc.

* **Stochastic effects** are malignant diseases and heritable effects for which the probability of an effect occurring, but not its severity, is a function of dose without threshold.

* **Lethal dose:** The whole-body dose from penetrating ionizing radiation resulting in 50% mortality in 30 days (assuming no medical treatment) is 2.5–4.5 Gy (250–450 rad)<fnref n="1" />, as measured internally on the body longitudinal center line. The surface dose varies due to variable body attenuation and may be a strong function of energy.

* **Cancer induction:**

The cancer induction probability is about 5% per Sv on average for the entire population [3].

* **Recommended effective dose limits:** The International Commission on Radiological Protection (ICRP) recommends a limit for radiation workers of 20 mSv effective dose per year averaged over 5 years, with the provision that the dose should not exceed 50 mSv in any single year [3]. The limit in the EU-countries

and Switzerland is 20 mSv per year, in the U.S. it is 50 mSv per year (5 rem per year). Many physics laboratories in the U.S. and elsewhere set lower limits. The effective dose limit for the general public is typically 1 mSv per year.

## 37.4 Prompt neutrons at accelerators

Neutrons dominate the radiation environment outside thick shielding (e.g., > 1 m of concrete) for high energy (> a few hundred MeV) electron and hadron accelerators. In addition, for accelerators with energies above about 10 GeV, muons contribute significantly at small angles with regard to the beam, even behind several meters of shielding. Another special case are synchrotron light sources where particular care has to be taken to shield the very intense low-energy photons extracted from the electron synchrotron into the experimental areas. Due to its importance at high energy accelerators this section focuses on prompt neutrons.

### 37.4.1 Electron accelerators

At electron accelerators, neutrons are generated via photonuclear reactions from bremsstrahlung photons. Neutron production takes place above a threshold value which varies from 10 to 19 MeV for light nuclei (with important exceptions, such as 2.23 MeV for deuterium and 1.67 MeV for beryllium) and from 4 to 6 MeV for heavy nuclei. It is commonly described by different mechanisms depending on the photon energy: the giant dipole resonance interactions (from threshold up to about 30 MeV, often the dominant process), the quasi-deuteron effect (between 30 MeV and a few hundred MeV), the delta resonance mechanism (between 200 MeV and a few GeV) and the vector meson dominance model at higher energies.

The giant dipole resonance reaction consists in a collective excitation of the nucleus, in which neutrons and protons oscillate in the direction of the photon electric field. The oscillation is damped by friction in a few cycles, with the photon energy being transferred to the nucleus in a process similar to evaporation. Nucleons emitted in the dipolar interaction have an anisotropic angular distribution, with a maximum at 90°, while those leaving the nucleus as a result of evaporation are emitted isotropically with a Maxwellian energy distribution described as [8]:

$$ \frac{dN}{dE_n} = \frac{E_n}{T^2} e^{-E_n/T}, \qquad (37.3) $$

where *T* is a nuclear 'temperature' (in units of MeV) characteristic of the particular target nucleus and its excitation energy. For heavy nuclei the 'temperature' generally lies in the range of *T* = 0.5–1.0 MeV. Neutron yields from semi-infinite targets per kW of electron beam power are plotted in Fig. 37.2 as a function of the electron beam energy [8].

While for thick targets neutron production is mainly due to photonuclear interactions, for thin targets (thickness of fractions of the radiation length) electronnuclear interactions are the dominating process.

Typical neutron energy spectra outside of concrete (80 cm thick, 2.35 g/cm³) and iron (40 cm thick) shields are shown in Fig. 37.3. In order to compare these spectra to those caused by proton beams (see below) the spectra are scaled by a factor of 100, which roughly corresponds to the difference in the high energy hadronic cross sections for photons and hadrons (e.g., the fine structure constant). The shape of these spectra are generally characterized by a low-energy peak at around 1 MeV (evaporation neutrons) and a high-energy shoulder at around 70–80 MeV. In case of concrete shielding, the spectrum also shows a pronounced peak at thermal neutron energies.

### 37.4.2 Proton accelerators

At proton accelerators, neutron yields emitted per incident proton by different target materials are roughly independent of proton energy between 20 MeV and 1 GeV, and are given by the ratio $C : Al : Cu-Fe : Sn : Ta-Pb = 0.3 : 0.6 : 1.0 : 1.5 : 1.7$ [11]. Above about 1 GeV, the neutron yield is proportional to $E^m$, where $0.80 \le m \le 0.85$ [12].

Typical neutron energy spectra outside of concrete and iron shielding are shown in Fig. 37.3. Here, the radiation fields are

<footnote n="1">¹RBE-weighted when necessary</footnote>
<figure><img src="image_252.png" /><figcaption>**Figure 37.2:** Neutron yields from semi-infinite targets per kW of electron beam power, as a function of the electron beam energy, disregarding target self-shielding [8].</figcaption></figure>

<figure><img src="image_253.png" /><figcaption>**Figure 37.3:** Neutron energy spectra calculated with the FLUKA code [9,10] from 25 GeV proton and electron beams on a thick copper target. Spectra are evaluated at 90° to the beam direction behind 80 cm of concrete or 40 cm of iron. All spectra are normalized per beam particle. For better visualization, spectra for electron beam are multiplied by a factor of 100.</figcaption></figure>

<figure><img src="image_254.png" /><figcaption>**Figure 37.4:** The variation of the attenuation length for mono-energetic neutrons in concrete as a function of neutron energy [11].</figcaption></figure>

caused by a 25 GeV proton beam interacting with a thick copper target. The comparison of these spectra with those for an electron beam of the same energy reflects the difference in the hadronic cross sections between photons and hadrons above a few 100 MeV. Differences are increasing towards lower energies because of different interaction mechanisms. Furthermore, the slight shift in energy above about 100 MeV follows from the fact that the energies of the interacting photons are lower than 25 GeV. Apart from this the shapes of the two spectra are similar.

The neutron-attenuation length is shown in Fig. 37.4 for concrete and mono-energetic broad-beam conditions. It reaches an asymptotic value of about 117 g/cm² above 200 MeV. As the cascade through thick shielding is carried by particles with energies between about 100 MeV and 300 MeV (in this energy range non-elastic cross sections are at minimum and are dominated by quasi-elastic processes leading to low attenuation) this value is equal to the equilibrium attenuation length for particles emitted at 90° in concrete.

## 37.5 Photon sources

The dose equivalent rate in tissue (in mSv/h) from a gamma point source emitting one photon of energy *E* (in MeV) per second at a distance of 1 m is $4.6 \times 10^{-9} \mu_{en}/\rho E$, where $\mu_{en}/\rho$ is the mass energy absorption coefficient. The latter has a value of 0.029 ± 0.004 cm²/g for photons in tissue over an energy range between 60 keV and 2 MeV (see Ref. [13] for tabulated values).

Similarly, the dose equivalent rate in tissue (in mSv/h) at the surface of a semi-infinite slab of uniformly activated material containing 1 Bq/g of a gamma emitter of energy *E* (in MeV) is $2.9 \times 10^{-4} R_{\mu} E$, where $R_{\mu}$ is the ratio of the mass energy absorption coefficients of the photons in tissue and in the material.

## 37.6 Accelerator-induced radioactivity

Typical medium- and long-lived activation products in metallic components of accelerators are $^{22}$Na, $^{46}$Sc, $^{48}$V, $^{51}$Cr, $^{54}$Mn, $^{55}$Fe, $^{59}$Fe, $^{56}$Co, $^{57}$Co, $^{58}$Co, $^{60}$Co, $^{63}$Ni and $^{65}$Zn. Gamma-emitting nuclides dominate doses from external irradiation at longer decay times (more than one day) while at short decay times $\beta^{+}$ emitters are also important (through photons produced by $\beta^{+}$ annihilation). Due to their short range, $\beta^{-}$ emitters are relevant, for example, only for dose to the skin and eyes or for doses due to inhalation or ingestion. Fig. 37.5 and Fig. 37.6 illustrate the contributions of gamma and $\beta^{+}$ emitters to the total dose rate at 12.4 cm distance to a copper sample [14]. The sample was activated by the stray radiation field created by a 120 GeV mixed hadron beam dumped in a copper target during about 8 hours at intensities between $10^{7} - 10^{8}$ hadrons per second. Here, the contributions by individual nuclides were calculated analytically assuming exponential decay of the nuclide inventory present at 20 minute cooling time. The total dose rate, however, also includes decay chains and as such also contributions from daughter nuclides.

Typically, dose rates at a certain decay time are mainly determined by radionuclides having a half-life of the order of the decay time. Extended irradiation periods might be an exception to this general rule as in this case the activity of long-lived nuclides can build up sufficiently so that it dominates that one of short-lived even at short cooling times.

Activation in concrete is dominated by $^{24}$Na (short decay time) and $^{22}$Na (long decay time). Both nuclides can be produced either by low-energy neutron reactions on the sodium-component in the concrete or by spallation reactions on silicon, calcium and other constituents such as aluminum. At long decay times nuclides of radiological interest in activated concrete can also be $^{60}$Co, $^{152}$Eu, $^{154}$Eu and $^{134}$Cs, all of which produced by (n,γ)-reactions with traces of natural cobalt, europium and cesium. Thus, such trace elements might be important even if their content in concrete is only a few parts per million or less by weight.

The explicit simulation of radionuclide production with general-purpose Monte Carlo codes has become the most commonly applied method to calculate induced radioactivity and its radiological consequences [14] (see also Sec. 37.8). They are complemented by analytical codes based on folding particle fluence spectra with nuclide production cross sections. ActiWiz [15,16] is an example
<figure><img src="image_255.png" /><figcaption>**Figure 37.5:** Contribution of individual gamma-emitting nuclides to the total dose rate at 12.4 cm distance to an activated copper sample [14].</figcaption></figure>

<figure><img src="image_256.png" /><figcaption>**Figure 37.6:** Contribution of individual positron-emitting nuclides to the total dose rate at 12.4 cm distance to an activated copper sample [14].</figcaption></figure>

of such a code targeting the domain of radiological characterization and material optimization. It allows for calculating nuclide inventories by convolution of fluence spectra with nuclide production data for 85 chemical elements and arbitrary compounds from threshold to an energy of 100 TeV.

## 37.7 Radiation protection instrumentation

The capacity to distinguish and measure the high-LET (mostly neutrons) and the low-LET components (photons, electrons, muons) of the radiation field at workplaces is of primary importance to evaluate the exposure of personnel. At proton accelerators the prompt dose equivalent outside a shield is mainly due to neutrons, with some contribution from photons and, to a minor extent, charged particles. At high-energy electron accelerators the dominant stray radiation during operation consists also of neutrons, because the shielding is normally thick enough to absorb most of the bremsstrahlung photons. Most of the personnel exposure at accelerator facilities is often received during maintenance interventions, and is due to gamma/beta radiation coming from

residual radioactivity in accelerator components.

Radiation detectors used both for radiation surveys and area monitoring are normally calibrated in ambient dose equivalent $H^*(10)$.

### 37.7.1 Neutron detectors

* **Rem counters:** A rem counter [17] is a portable detector consisting of a thermal neutron counter embedded in a polyethylene moderator, with a response function that approximately follows the curve of the conversion coefficients from neutron fluence to $H^*(10)$ over a wide energy range. Conventional rem counters provide a response to neutrons up to approximately 10-15 MeV, extended-range units are heavier as they include a converter consisting of a material with high mass number to correctly measure $H^*(10)$ up to several hundred MeV [18].

* **Bonner Sphere Spectrometer (BSS):** A BSS [19] consists of a thermal neutron detector at the centre of moderating spheres of different diameters made of polyethylene (PE) or a combination of PE and a high-A material to enhance its response to high energy neutrons (similar to rem counters). Each sphere has a different response function versus neutron energy, and the neutron energy, at which the sensitivity peaks, increases with sphere diameter. The energy resolution of the system is rather low but satisfactory for radiation protection purposes. The neutron spectrum is obtained by unfolding the experimental counts of the BSS with its response matrix by a computer code that is often based on an iterative algorithm such as GRAVEL [20] and MAXED [21]. BSS exist in active (using $^3$He or $BF_3$ proportional counters or $^6$LiI scintillators) and passive versions (using CR-39 track detectors or LiF), for use *e.g.* in strongly pulsed fields. With $^3$He counters the discrimination with respect to gamma rays and noise is excellent.

* **Bubble detectors:** A bubble detector [22] is a dosimeter based on a super-heated emulsion (super-heated droplets suspended in a gel) contained in a vial and acting as a continuously sensitive, miniature bubble chamber. The total number of bubbles evolved from the radiation-induced nucleation of drops gives an integrated measure of the total neutron exposure. Various techniques exist to record and count the bubbles, *e.g.*, visual inspection, automated reading with video cameras or acoustic counting. Bubble detectors are insensitive to low-LET radiation. Super-heated emulsions are used as personal, area and environmental dosimeters, as well as neutron spectrometers.

* **Track etched detectors:** Track etched detectors (TEDs) [23] are based on the preferential dissolution of suitable, mostly insulator, materials along the damage trails of charged particles of sufficiently high-energy deposition density. The detectors are effectively not sensitive to radiation which deposits the energy through the interactions of particles with low LET. These dosimeters are generally able to determine neutron ambient dose equivalent down to around 100 µSv. They are used both as personal dosimeters and for area monitoring, *e.g.*, in BSS.

### 37.7.2 Photon detectors [24]

* **Geiger Müller (GM) counters:** GM counters are low cost devices simple to operate. They work in pulse mode and since they only count radiation-induced events, any spectrometric information is lost. In general they are calibrated in terms of air kerma, for instance in a $^{60}$Co field. The response of GM counters to photons is constant within 15% for energies up to 2 MeV and shows considerable energy dependence above.

* **Ionization chambers:** Ionization chambers are gas-filled detectors used both as hand-held instruments (*e.g.*, for radiation surveys) and environmental monitors. They are normally operated in current mode although pulse-mode operation is also possible. They possess a relatively flat response to a wide range of X- and gamma ray energies (typically from 10 keV to several MeV), can measure radiation over a wide intensity range and are capable of discriminating between the beta and gamma components of a radiation field (by use of, *e.g.*, a beta window). Pressurized ion chambers (filled, *e.g.*, with Ar or H gas to several tens of bars) are used for environmental and stray radiation monitoring applications. They have good sensitivity to neutrons and charged hadrons in addition to low LET radiation (gammas and muons), with the response function to the former being strongly non-linear
with energy.

* **Proportional counters:** Proportional counters are another type of gas-filled detector, operated in pulse mode. They are normally of cylindrical geometry, with a thin central wire (the anode) on the axis of a hollow tube (the cathode). The initial electron charge created by the interaction of the radiation with the gas is amplified by the multiplication process (Townsend avalanche) produced by the strong electric field in a narrow volume around the anode. The generated pulses are still proportional to the original number of ion pairs but are much larger than in an ionization chamber, yielding a much improved signal-to-noise ratio. Proportional counters are also used for neutron detection: the rem counter and BSS discussed in the previous sections typically employ ³He or BF₃ filled tubes with a large sensitivity to thermal neutrons.

* **Scintillators:** Scintillation-based detectors are used in radiation protection as hand-held probes and in fixed installations, e.g., portal monitors. A scintillation detector or counter is obtained coupling a scintillator to an electronic light sensor such as a photomultiplier tube (PMT), a photodiode or a silicon photomultiplier (SiPM). There is a wide range of scintillating materials, inorganic (such as CsI and BGO), organic or plastic; they find application in both photon dosimetry and spectrometry.

### 37.7.3 Operation in pulsed radiation fields

There are many practical situations with particle accelerators used for scientific, industrial and medical applications where the time structure of the stray radiation limits the use of active monitors or requires specifically designed electronics. Pulsed neutron and gamma fields may be present because of beam losses at, e.g., targets, collimators and beam dumps. The time duration of a single burst can range from a few ns to about 1 ms with a typical repetition rate in the range 0.1–100 Hz. Conventional detectors generally suffer from dead time effects and have strong limitations in the measurements of pulsed fields. Severe under response has been observed, e.g., in commercial rem counters, with tremendous underestimation of the ambient dose equivalent, $H^*(10)$, up to three orders of magnitude [25]. The common techniques used to correct the response of radiation detectors which include dead-time corrections operate properly in a steady-state radiation field, whereas it is much more difficult to cope with dead time losses in a pulsed radiation field of unknown time structure and burst dose. Generally speaking the monitoring instrumentation must be chosen on the basis of the knowledge of the radiation field. Detectors with specifically designed electronics must be employed in pulsed field conditions in place of conventional rem counters for neutrons, such as the LUPIN detector [26]. If real-time monitoring is not required, passive detectors or dosimeters such as TEDs mentioned in Sec. 37.7.1 or LiF mentioned in Sec. 37.7.5 can be employed, as they are insensitive to the time structure of the radiation.

### 37.7.4 Operation in presence of magnetic fields

There may be circumstances in which a radiation survey needs to be performed in the presence of a magnetic field. This may be the case for example for measurements on activated detector components in the underground experimental areas of CERN LHC experiments without switching off the intense magnetic field of the detectors; or measurements of the residual radioactivity of permanent magnets sometimes used in particle accelerators. In hospitals, one of most advanced diagnostic modalities is PET-MRI, in which PET imaging uses radioactive substances (radio-pharmaceuticals) in proximity of the strong magnetic field of the MRI scanner. Most of the commercial survey meters do not operate correctly or do not function at all in the presence of even comparatively weak magnetic fields below 0.3 T, and often already at 0.1 T. This is either because of the presence of components in their electronics sensitive to magnetic field (e.g. inductive elements), ferromagnetic materials in general (such as shells, screws and frames), battery drain/heating or LCD screen failure. An instrument called B-RAD has recently been developed [27] that overcomes the above limitations and is now commercially available. It is a hand-held instrument for gamma dose rate and gamma spectrometry measurements, consisting of a central unit (housing a microcontroller, magnetic field insensitive electronic

circuitry, the battery and two displays) and a probe connected through a cable, certified to operate up to 3 T fields. The probe uses a LaBr₃(Ce) crystal coupled to an array of silicon photomultipliers. The excellent scintillation properties and the high photon resolution of the scintillator (3.3% FWHM at 662 keV) make the device capable of operating in the energy range from 30 keV to 2 MeV photons with a very fast response.

### 37.7.5 Personal dosimeters

Personal dosimeters, calibrated in $H_p(10)$, are worn by persons exposed to ionizing radiation for professional reasons to record the dose received. They are typically passive detectors, either film, track etched detectors, $^6$Li/$^7$Li-based dosimeters (e.g. LiF), optically stimulated luminescence (OSL) or radiophotoluminescence detectors (RPL) but semi-active dosimeters using miniaturized ion-chambers also exist, like the Direct Ion Storage (DIS) dosimeters in use at CERN.

Electronic personal dosimeters are small active units for on-line monitoring of individual exposure, designed to be worn on the body. They can give an alarm on both the integral dose received or dose rate once a pre-set threshold is exceeded.

## 37.8 Monte Carlo codes for radiation protection studies

The use of general-purpose particle interaction and transport Monte Carlo codes is often the most accurate and efficient choice for assessing radiation protection quantities at accelerators. Due to the vast spread of such codes to all areas of particle physics and the associated extensive benchmarking with experimental data, the modeling has reached an unprecedented accuracy. Furthermore, most codes allow the user to simulate all aspects of a high energy particle cascade in one and the same run: from the first interaction of a TeV nucleus over the transport and re-interactions (hadronic and electromagnetic) of the produced secondaries, to detailed nuclear fragmentation, the calculation of radioactive decays and even of the electromagnetic shower caused by the radiation from such decays. A brief account of the codes most widely used for radiation protection studies at high energy accelerators is given in the following.

* **FLUKA [9, 10, 28, 29]:** FLUKA is a general-purpose particle interaction and transport code. It comprises all features needed for radiation protection, such as detailed hadronic and nuclear interaction models up to 10 PeV, full coupling between hadronic and electromagnetic processes and numerous variance reduction options. The latter include weight windows, region importance biasing, and leading particle, interaction, and decay length biasing (among others). The capabilities of FLUKA are unique for studies of induced radioactivity, especially with regard to nuclide production, decay, and transport of residual radiation. In particular, particle cascades by prompt and residual radiation are simulated in parallel based on the microscopic models for nuclide production and a solution of the Bateman equations for activity build-up and decay.

* **GEANT4 [30–32]:** GEANT4 is an object-oriented toolkit consisting of a kernel that provides the framework for particle transport, including tracking, geometry description, material specifications, management of events and interfaces to external graphics systems. The kernel also provides interfaces to physics processes. It allows the user to freely select the physics models that best serve the particular application needs. Implementations of interaction models exist over an extended range of energies, from optical photons and thermal neutrons to high-energy interactions required for the simulation of accelerator and cosmic ray experiments. To facilitate the use of variance reduction techniques, general-purpose biasing methods such as importance biasing, weight windows, and a weight cut-off method have been introduced directly into the toolkit. Other variance reduction methods, such as leading particle biasing for hadronic processes, come with the respective physics packages.

* **MARS15 [33,34]:** The MARS15 code system is a set of Monte Carlo programs for the simulation of hadronic and electromagnetic cascades. It covers a wide energy range: 1 keV to 100 TeV for muons, charged hadrons, heavy ions and electromagnetic showers; and 0.00215 eV to 100 TeV for neutrons. Hadron-nucleus inter-
actions as well as practically all other strong, weak and electro-
magnetic interactions in the entire energy range can be simulated
either inclusively or exclusively. MARS15 uses ENDFB-VII nu-
clear data to handle interactions of neutrons with energies below
14 MeV. Several variance reduction techniques, such as weight
windows, particle splitting, and Russian roulette, are available.
A tagging module allows tagging the origin of a given signal for
source term or sensitivity analyses. The geometry module al-
lows either a basic solid body representation option or a ROOT-
based powerful engine. Further features of MARS15 include a
MAD-MARS merge for a convenient creation of accelerator mod-
els and multi-turn tracking and cascade simulation in accelerator
and beamline lattices.

* MCNP6 [35,36]: MCNP6 is the latest version of the Monte Carlo N-Particle transport (MCNP) family of neutron interaction and transport codes and, therefore, features one of the most comprehensive and detailed descriptions of the related physical processes. It transports 37 different particle types, including ions and electromagnetic particles. The neutron interaction and transport modules use standard evaluated data libraries mixed with physics models where such libraries are not available. The transport is continuous in energy. MCNP6 contains one of the most powerful implementations of variance reduction techniques. Spherical mesh weight windows can be created by a generator in order to focus the simulation time on certain spatial regions of interest. In addition, a more generalized phase space biasing is also possible through energy- and time-dependent weight windows. Other biasing options include pulse-height tallies with variance reduction and criticality source convergence acceleration.

* PHITS [37,38]: The Particle and Heavy-Ion Transport code System PHITS was among the first general-purpose codes to simulate the transport and interactions of heavy ions in a wide energy range, from 10 MeV/nucleon to 100 GeV/nucleon. It is based on the high-energy hadron transport code NMTC/JAM that was extended to heavy ions. The transport of low-energy neutrons employs cross sections from evaluated nuclear data libraries such as ENDF and JENDL below 20 MeV. Electromagnetic interactions are simulated based on the ITS code in the energy range between 1 keV and 100 MeV for electrons and positrons and between 1 keV and 100 GeV for photons. Several variance reduction techniques, including weight windows and region importance biasing, are available.

References

[1] International Commission on Radiation Units and Measure-
ments, *Fundamental Quantities and Units for Ionizing Ra-
diation*, ICRU Report 60 (1998).

[2] ICRP, 2010. *Conversion Coefficients for Radiological Pro-
tection Quantities for External Radiation Exposures*, ICRP
Publication 116, Annals of the ICRP 40(2-5).

[3] ICRP Publication 103, *The 2007 Recommendations of the International Commission on Radiological Protection*, Annals of the ICRP, Elsevier (2007).

[4] M. Pelliccioni, Radiation Protection Dosimetry **88**, 279 (2000).

[5] E. Pochin, *Nuclear Radiation: Risks and Benefits*, Clarendon Press, Oxford, 1983.

[6] United Nations, *Report of the United Nations Scientific Committee on the Effect of Atomic Radiation*, General Assembly, Official Records A/63/46 (2008).

[7] G. Cinelli et al., *European Atlas of Natural Radiation*, G. Cinelli, M. De Cort, and T. Tollefsen, editor(s), Publications Office of the European Union, Luxembourg, 2019, ISBN 978-92-76-08259-0, JRC116795.

[8] W. P. Swanson, *Radiological Safety Aspects of the Operation of Electron Linear Accelerators*, IAEA Technical Reports Series No. 188 (1979) (1979).

[9] A. Ferrari, et al., FLUKA: A Multi-particle Transport Code, CERN-2005-010 (2005), SLAC-R-773, INFN-TC-05-11, http://www.fluka.org.

[10] G. Battistoni et al., Annals of Nuclear Energy **82**, 10 (2015).

[11] R.H. Thomas and G.R. Stevenson, *Radiological Safety Aspects of the Operation of Proton Accelerators*, SSC-N-354 (1986), IAEA Technical Report Series No. 283 (1988).

[12] T. Gabriel et al., Nucl. Instrum. Meth. A **338**, 336 (1994).

[13] http://physics.nist.gov/PhysRefData/XrayMassCoef/cover.html.

[14] S. Roesler, et al., "Simulation of Remanent Dose Rates and Benchmark Measurements at the CERN-EU High Energy Reference Field Facility," in *Proceedings of the Sixth International Meeting on Nuclear Applications of Accelerator Technology*, San Diego, CA, 1-5 June 2003, 655-662 (2003).

[15] C. Theis and H. Vincke, “The use of ActiWiz in opera-
tional radiation protection,” in *Proceedings of the Twelfth
Meeting of Task-Force on Shielding Aspects of Accelerators,
Targets and Irradiation Facilities of Accelerator Technology,
SATIF12* FNAL, 28-30 April 2014, Nuclear Science Report
NEA/NSC/R 3, (2015).

[16] C. Theis and H. Vincke, “ActiWiz3 an overview of the concepts, architecture and new features,” CERN Technical Note CERN-RP-2016-117-REPORTS-TN (2016), http://actiwiz.web.cern.ch .

[17] I.O. Andersson and J. Braun, “A neutron rem counter with uniform sensitivity from 0.025 eV to 10 MeV” in Proceedings of the IAEA Symposium on Neutron dosimetry, IAEA, Vienna, Vol. II, 87–95, (1963).

[18] C. Birattari, et al., Radiation Protection Dosimetry **76**, 135 (1998).

[19] R.L. Bramblett, R.I. Ewing and T.W. Bonner, Nucl. Instrum. Methods **9**, 1 (1960).

[20] M. Matzke, PTB, Braunschweig PTBN-19, (1994).

[21] M. Reginatto and P. Goldhagen, Health Physics **77**, 579 (1999).

[22] F. d'Errico, Nucl. Instrum. Methods **B184**, 229 (2001).

[23] K. Becker, Dosimetric applications of track etching, in *Topics in Radiation Dosimetry* Ed. F.H. Attix, Academic Press, London, 79-143, (1972).

[24] G.F. Knoll, *Radiation detection and measurements*, Wiley (2010).

[25] M. Caresana et al., Nucl. Instrum. Meth. A **737**, 203 (2014).

[26] M. Caresana, et al., A new version of the LUPIN detector: Improvements and latest experimental verification, Review of Scientific Instruments **85**, 065102 (2014).

[27] D. Celeste, A. Curioni, A. Fazzi, M. Silari and V. Varoli, B-RAD: a radiation survey meter for operation in intense magnetic fields Journal of Instrumentation **14**, T05007 (2019).

[28] C. Ahdida, et al., New Capabilities of the FLUKA Multi-Purpose Code, Frontiers in Physics **9**, 788253 (2022) http://fluka.cern .

[29] T. T. Böhlen et al., Nucl. Data Sheets **120**, 211 (2014).

[30] S. Agostinelli et al. (GEANT4), Nucl. Instrum. Meth. A **506**, 250 (2003).

[31] J. Allison et al., IEEE Trans. Nucl. Sci. **53**, 270 (2006).

[32] J. Allison et al., Nucl. Instrum. Meth. A **835**, 186 (2016), http://geant4.cern.ch.

[33] N.V. Mokhov and C.C. James, The MARS Code System User's Guide, Fermilab-FN-1058-APC (2018), https://mars.fnal.gov.

[34] N. Mokhov et al., Prog. Nucl. Sci. Tech. **4**, 496 (2014), [arXiv:1409.0033].

[35] J.T. Goorley, et al., Nuclear Technology **180**, 298 (2012).

[36] J.A. Kulesza et al., MCNP Code Version 6.3.0 Theory and User Manual, Los Alamos National Laboratory report, LA-UR-22-30006 (2022), https://mcnp.lanl.gov.

[37] T. Sato, et al., Journal of Nuclear Science and Technology **55**, 684 (2018).

[38] T. Sato, et al., PHITS Particle and Heavy Ion Transport code System, Version 3.31 (2023), https://phits.jaea.go.jp.
Mathematical Tools or Statistics, Monte Carlo, Group Theory

<table><tr><td>39. Probability</td><td>689</td></tr><tr><td>40. Statistics (rev.)</td><td>694</td></tr><tr><td>41. Machine Learning</td><td>711</td></tr><tr><td>42. Monte Carlo techniques</td><td>748</td></tr><tr><td>43. Monte Carlo event generators (rev.)</td><td>752</td></tr><tr><td>44. Monte Carlo neutrino event generators</td><td>764</td></tr><tr><td>45. Monte Carlo particle numbering scheme (rev.)</td><td>768</td></tr><tr><td>46. Clebsch-Gordan coefficients, spherical harmonics, and d functions</td><td>772</td></tr><tr><td>47. SU(3) isoscalar factors and representation matrices</td><td>773</td></tr><tr><td>48. SU(n) multiplets and Young diagrams</td><td>774</td></tr></table>
<figure><img src="image_257.png" /></figure>
# 39. Probability

Revised August 2021 by G. Cowan (RHUL).

## 39.1 General

Further discussion of probability can be found, e.g., in Refs. [1–8].

An abstract definition of probability can be given by considering a set $S$, called the sample space, and possible subsets $A, B, \dots$, the interpretation of which is left open. The probability $P$ is a real-valued function defined by the following axioms due to Kolmogorov [9]:

1. For every subset $A$ in $S$, $P(A) \ge 0$;

2. For disjoint subsets (i.e., $A \cap B = \emptyset$), $P(A \cup B) = P(A) + P(B)$;

3. $P(S) = 1$.

In addition, one defines the conditional probability $P(A|B)$ (read as $P$ of $A$ given $B$) as

$$P(A|B) = \frac{P(A \cap B)}{P(B)} \qquad (39.1)$$

From this definition and using the fact that $A \cap B$ and $B \cap A$ are the same, one obtains Bayes' theorem,

$$P(A|B) = \frac{P(B|A)P(A)}{P(B)} \qquad (39.2)$$

From the three axioms of probability and the definition of conditional probability, one obtains the law of total probability,

$$P(B) = \sum_{i} P(B|A_{i})P(A_{i}) \qquad (39.3)$$

for any subset $B$ and for disjoint $A_i$ with $\cup_i A_i = S$. This can be combined with Bayes' theorem (Eq. (39.2)) to give

$$P(A|B) = \frac{P(B|A)P(A)}{\sum_{i} P(B|A_{i})P(A_{i})} \qquad (39.4)$$

where the subset $A$ could, for example, be one of the $A_i$.

The most commonly used interpretation of the elements of the sample space are outcomes of a repeatable experiment. The probability $P(A)$ is assigned a value equal to the limiting frequency of occurrence of $A$. This interpretation forms the basis of frequentist statistics.

The elements of the sample space might also be interpreted as hypotheses, i.e., statements that are either true or false, such as ‘The mass of the W boson lies between 80.3 and 80.5 GeV’. Upon repetition of a measurement, however, such statements are either always true or always false, i.e., the corresponding probabilities in the frequentist interpretation are either 0 or 1. Using subjective probability, however, $P(A)$ is interpreted as the degree of belief that the hypothesis $A$ is true. Subjective probability is used in Bayesian (as opposed to frequentist) statistics. Bayes' theorem can be written

$$P(\text{theory}|\text{data}) \propto P(\text{data}|\text{theory})P(\text{theory}) \qquad (39.5)$$

where ‘theory’ represents some hypothesis and ‘data’ is the outcome of the experiment. Here $P(\text{theory})$ is the prior probability for the theory, which reflects the experimenter’s degree of belief before carrying out the measurement, and $P(\text{data}|\text{theory})$ is the probability to have gotten the data actually obtained, given the theory, which is also called the likelihood.

Bayesian statistics provides no fundamental rule for obtaining the prior probability, which may depend on previous measurements, theoretical prejudices, etc. Once this has been specified, however, Eq. (39.5) tells how the probability for the theory must be modified in the light of the new data to give the posterior probability, $P(\text{theory}|\text{data})$. As Eq. (39.5) is stated as a proportionality, the probability must be normalized by summing (or integrating) over all possible hypotheses.

## 39.2 Random variables

A random variable is a numerical characteristic assigned to an element of the sample space. In the frequency interpretation of probability, it corresponds to an outcome of a repeatable experiment. Let $x$ be a possible outcome of an observation. If $x$ can take on any value from a continuous range, we write $f(x; \theta)dx$ as the probability that the measurement's outcome lies between $x$ and $x + dx$. The function $f(x; \theta)$ is called the probability density function (p.d.f.), which may depend on one or more parameters $\theta$. If $x$ can take on only discrete values (e.g., the non-negative integers), then we use $f(x; \theta)$ to denote the probability to find the value $x$. In the following the term p.d.f. is often taken to cover both the continuous and discrete cases, although technically the term density should only be used in the continuous case.

The p.d.f. is always normalized to unity. Both $x$ and $\theta$ may have multiple components and are then often written as vectors. If $\theta$ is unknown, we may wish to estimate its value from a given set of measurements of $x$; this is a central topic of statistics (see Sec. 40).

The cumulative distribution function $F(a)$ is the probability that $x \le a$:

$$F(a) = \int_{-\infty}^{a} f(x) dx \qquad (39.6)$$

Here and below, if $x$ is discrete-valued, the integral is replaced by a sum. The endpoint $a$ is expressly included in the integral or sum. Then $0 \le F(x) \le 1$, $F(x)$ is nondecreasing, and $P(a < x \le b) = F(b) - F(a)$. If $x$ is discrete, $F(x)$ is flat except at allowed values of $x$, where it has discontinuous jumps equal to $f(x)$.

Any function of random variables is itself a random variable, with (in general) a different p.d.f. The expectation value of any function $u(x)$ is

$$E[u(x)] = \int_{-\infty}^{\infty} u(x) f(x) dx \qquad (39.7)$$

assuming the integral is finite. The expectation value is linear, i.e., for any two functions $u$ and $v$ of $x$ and constants $c_1$ and $c_2$, $E[c_1u + c_2v] = c_1E[u] + c_2E[v]$.

The $n^{th}$ moment of a random variable $x$ is

$$\alpha_n \equiv E[x^n] = \int_{-\infty}^{\infty} x^n f(x) dx \qquad (39.8a)$$

and the $n^{th}$ central moment of $x$ (or moment about the mean, $\alpha_1$) is

$$m_n \equiv E[(x - \alpha_1)^n] = \int_{-\infty}^{\infty} (x - \alpha_1)^n f(x) dx \qquad (39.8b)$$

The most commonly used moments are the mean $\mu$ and variance $\sigma^2$:

$$\mu \equiv \alpha_1 \qquad (39.9a)$$

$$\sigma^2 \equiv V[x] \equiv m_2 = \alpha_2 - \mu^2 \qquad (39.9b)$$

The mean is the location of the “center of mass” of the p.d.f., and the variance is a measure of the square of its width. Note that $V[cx + k] = c^2V[x]$. It is often convenient to use the standard deviation of $x$, $\sigma$, defined as the square root of the variance.

Any odd moment about the mean is a measure of the skewness of the p.d.f. The simplest of these is the dimensionless coefficient of skewness $\gamma_1 = m_3/\sigma^3$.

The fourth central moment $m_4$ provides a convenient measure of the tails of a distribution. For the Gaussian distribution (see Sec. 39.4), one has $m_4 = 3\sigma^4$. The kurtosis is defined as $\gamma_2 = m_4/\sigma^4 - 3$, i.e., it is zero for a Gaussian, positive for a lep-tokurtic distribution with longer tails, and negative for a platykur-tic distribution with tails that die off more quickly than those of a Gaussian.

The quantile $x_\alpha$ is the value of the random variable $x$ at which the cumulative distribution is equal to $\alpha$. That is, the quantile
is the inverse of the cumulative distribution function, i.e., $x_α = F^{-1}(α)$. An important special case is the median, $x_{med}$, defined by $F(x_{med}) = 1/2$, i.e., half the probability lies above and half lies below $x_{med}$. (More rigorously, $x_{med}$ is a median if $P(x \ge x_{med}) \ge 1/2$ and $P(x \le x_{med}) \ge 1/2$. If only one value exists, it is called 'the median.')

Under a monotonic change of variable $x \to y(x)$, the quantiles of a distribution (and hence also the median) obey $y_α = y(x_α)$. In general the expectation value and mode (most probable value) of a distribution do not, however, transform in this way.

Let $x$ and $y$ be two random variables with a joint p.d.f. $f(x, y)$. The marginal p.d.f. of $x$ (the distribution of $x$ with $y$ unobserved) is

$$f_1(x) = \int_{-\infty}^{\infty} f(x, y) dy, \qquad (39.10)$$

and similarly for the marginal p.d.f. $f_2(y)$. The conditional p.d.f. of $y$ given fixed $x$ (with $f_1(x) \neq 0$) is defined by $f_3(y|x) = f(x, y)/f_1(x)$, and similarly $f_4(x|y) = f(x, y)/f_2(y)$. From these, we immediately obtain Bayes' theorem (see Eqs. (39.2) and (39.4)),

$$f_4(x|y) = \frac{f_3(y|x)f_1(x)}{f_2(y)} = \frac{f_3(y|x)f_1(x)}{\int f_3(y|x')f_1(x')dx'}. \quad (39.11)$$

The mean of $x$ is

$$\mu_x = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x f(x, y) dx dy = \int_{-\infty}^{\infty} x f_1(x) dx, \quad (39.12)$$

and similarly for $y$. The covariance of $x$ and $y$ is

$$\operatorname{cov}[x, y] = E[(x - \mu_x)(y - \mu_y)] = E[xy] - \mu_x \mu_y. \quad (39.13)$$

A dimensionless measure of the covariance of $x$ and $y$ is given by the correlation coefficient,

$$\rho_{xy} = \operatorname{cov}[x, y] / \sigma_x \sigma_y, \quad (39.14)$$

where $\sigma_x$ and $\sigma_y$ are the standard deviations of $x$ and $y$. It can be shown that $-1 \le \rho_{xy} \le 1$.

Two random variables $x$ and $y$ are independent if and only if

$$f(x, y) = f_1(x)f_2(y). \quad (39.15)$$

If $x$ and $y$ are independent, then $\rho_{xy} = 0$; the converse is not necessarily true. If $x$ and $y$ are independent, $E[u(x)v(y)] = E[u(x)]E[v(y)]$, and $V[x+y] = V[x] + V[y]$; otherwise, $V[x+y] = V[x] + V[y] + 2\operatorname{cov}[x, y]$, and $E[uv]$ does not necessarily factorize.

Consider a set of $n$ continuous random variables $x = (x_1, \dots, x_n)$ with joint p.d.f. $f(x)$, and a set of $n$ new variables $y = (y_1, \dots, y_n)$, related to $x$ by means of a function $y(x)$ that is one-to-one, i.e., the inverse $x(y)$ exists. The joint p.d.f. for $y$ is given by

$$g(y) = f(x(y))|J|, \quad (39.16)$$

where $|J|$ is the absolute value of the determinant of the square matrix $J_{ij} = \partial x_i / \partial y_j$ (the Jacobian determinant). If the transformation from $x$ to $y$ is not one-to-one, the $x$-space must be broken into regions where the function $y(x)$ can be inverted, and the contributions to $g(y)$ from each region summed.

Given a set of functions $y = (y_1, \dots, y_m)$ with $m < n$, one can construct $n-m$ additional independent functions, apply the procedure above, then integrate the resulting $g(y)$ over the unwanted $y_i$ to find the marginal distribution of those of interest.

For a one-to-one transformation of discrete random variables, the probability is obtained by simple substitution; no Jacobian is necessary because in this case $f$ is a probability rather than a probability density. If the transformation is not one-to-one, then one must sum the probabilities for all values of the original variable that contribute to a given value of the transformed variable. If $f$ depends on a set of parameters $\theta$, a change to a different parameter set $\eta(\theta)$ is made by simple substitution; no Jacobian is used.

### 39.2.1 Propagation of errors

Consider $n$ random variables $x = (x_1, \dots, x_n)$ and $m$ functions $y(x) = (y_1(x), \dots, y_m(x))$. Suppose here that the mean values $\mu = (\mu_1, \dots, \mu_n) = E[x]$ are known, although in practice they will only be estimated, and suppose we also know or have estimated the covariance matrix $V_{ij} = \operatorname{cov}[x_i, x_j]$. The goal of error propagation is to determine the covariance matrix for the functions, $U_{ij} = \operatorname{cov}[y_i, y_j]$. In particular, the diagonal elements $U_{ii} = V[y_i]$ give the variances. The new covariance matrix can be found by expanding the functions $y(x)$ about the means $\mu$ to first order in a Taylor series. Using this one finds

$$U_{ij} \approx \sum_{k,l} \left. \frac{\partial y_i}{\partial x_k} \frac{\partial y_j}{\partial x_l} \right|_{\mu} V_{kl}. \quad (39.17)$$

This can be written in matrix notation as $U \approx AVA^T$ where the matrix of derivatives $A$ is

$$A_{ij} = \left. \frac{\partial y_i}{\partial x_j} \right|_{\mu}, \quad (39.18)$$

and $A^T$ is its transpose. The approximation is exact if $y(x)$ is linear. If this is not the case, the approximation can break down if, for example, $y(x)$ is significantly nonlinear close to $\mu$ in a region of a size comparable to the standard deviations of $x$.

## 39.3 Characteristic functions

The characteristic function $\phi(u)$ associated with the p.d.f. $f(x)$ is essentially its Fourier transform, or the expectation value of $e^{iux}$:

$$\phi(u) = E[e^{iux}] = \int_{-\infty}^{\infty} e^{iux} f(x) dx. \quad (39.19)$$

Once $\phi(u)$ is specified, the p.d.f. $f(x)$ is uniquely determined and vice versa; knowing one is equivalent to the other. Characteristic functions are useful in deriving a number of important results about moments and sums of random variables.

It follows from Eqs. (39.8) and (39.19) that the $n^{th}$ moment of a random variable $x$ that follows $f(x)$ is given by

$$i^{-n} \left. \frac{d^n \phi}{du^n} \right|_{u=0} = \int_{-\infty}^{\infty} x^n f(x) dx = \alpha_n. \quad (39.20)$$

Thus it is often easy to calculate all the moments of a distribution defined by $\phi(u)$, even when $f(x)$ cannot be written down explicitly.

If the p.d.f.s $f_1(x)$ and $f_2(y)$ for independent random variables $x$ and $y$ have characteristic functions $\phi_1(u)$ and $\phi_2(u)$, then the characteristic function of the weighted sum $ax + by$ is $\phi_1(au)\phi_2(bu)$. The rules of addition for several important distributions (e.g., that the sum of two Gaussian distributed variables also follows a Gaussian distribution) easily follow from this observation.

Let the (partial) characteristic function corresponding to the conditional p.d.f. $f_2(x|z)$ be $\phi_2(u|z)$, and the p.d.f. of $z$ be $f_1(z)$. The characteristic function after integration over the conditional value is

$$\phi(u) = \int \phi_2(u|z) f_1(z) dz. \quad (39.21)$$

Suppose we can write $\phi_2$ in the form

$$\phi_2(u|z) = A(u)e^{ig(u)z}. \quad (39.22)$$

Then

$$\phi(u) = A(u)\phi_1(g(u)). \quad (39.23)$$

The cumulants (semi-invariants) $\kappa_n$ of a distribution with characteristic function $\phi(u)$ are defined by the relation

$$\phi(u) = \exp \left[ \sum_{n=1}^{\infty} \frac{\kappa_n}{n!} (iu)^n \right] = \exp \left( i\kappa_1 u - \frac{1}{2} \kappa_2 u^2 + \dots \right). \quad (39.24)$$
The values $\kappa_n$ are related to the moments $\alpha_n$ and $m_n$. The first few relations are

$$
\begin{align}
\kappa_1 &= \alpha_1 (= \mu, \text{ the mean}) \nonumber \\
\kappa_2 &= m_2 = \alpha_2 - \alpha_1^2 (= \sigma^2, \text{ the variance}) \tag{39.25} \\
\kappa_3 &= m_3 = \alpha_3 - 3\alpha_1\alpha_2 + 2\alpha_1^3. \nonumber
\end{align}
$$

39.4 Commonly used probability distributions

Table 39.1 gives a number of common probability density func-
tions and corresponding characteristic functions, means, and vari-
ances. Further information may be found in Refs. [1–8] [10],
and [11], which has particularly detailed tables. Monte Carlo tech-
niques for generating each of them may be found in our Sec. 42.4
and in Ref. [10]. We comment below on all except the trivial
uniform distribution.

39.4.1 Binomial and multinomial distributions

A random process with exactly two possible outcomes which
occur with fixed probabilities is called a *Bernoulli* process. If the
probability of obtaining a certain outcome (a “success”) in an in-
dividual trial is *p*, then the probability of obtaining exactly *r* suc-
cesses (*r* = 0, 1, 2, ..., *N*) in *N* independent trials, without regard
to the order of the successes and failures, is given by the binomial
distribution *f*(*r*; *N*, *p*) in Table 39.1. If *r* and *s* are binomially
distributed with parameters (*N*<sub>r</sub>, *p*) and (*N*<sub>s</sub>, *p*), then *t* = *r* + *s*
follows a binomial distribution with parameters (*N*<sub>r</sub> + *N*<sub>s</sub>, *p*). If
there are *m* possible outcomes for each trial having probabili-
ties *p*<sub>1</sub>, *p*<sub>2</sub>, ..., *p*<sub>m</sub>, then the joint probability to find *r*<sub>1</sub>, *r*<sub>2</sub>, ..., *r*<sub>m</sub>
of each outcome after a total of *N* independent trials is given by
the multinomial distribution as shown in Table 39.1. We can
regard outcome *i* as “success” and all the rest as “failure”, so in-
dividually, any of the *r*<sub>*i*</sub> follow a binomial distribution for *N* trials
and a success probability *p*<sub>*i*</sub>.

39.4.2 Poisson distribution

The Poisson distribution $f(n; \nu)$ gives the probability of finding exactly $n$ events in a given interval of $x$ (e.g., space or time) when the events occur independently of one another and of $x$ at an average rate of $\nu$ per the given interval. The variance $\sigma^2$ equals $\nu$. It is the limiting case $p \to 0$, $N \to \infty$, $Np = \nu$ of the binomial distribution. The Poisson distribution approaches the Gaussian distribution for large $\nu$.

For example, a large number of radioactive nuclei of a given
type will result in a certain number of decays in a fixed time
interval. If this interval is small compared to the mean lifetime,
then the probability for a given nucleus to decay is small, and
thus the number of decays in the time interval is well modeled as
a Poisson variable.

39.4.3 Normal or Gaussian distribution

The normal (or Gaussian) probability density function
$f(x; \mu, \sigma^2)$ given in Table 39.1 has mean $E[x] = \mu$ and variance
$V[x] = \sigma^2$. Comparison of the characteristic function $\phi(u)$ given
in Table 39.1 with Eq. (39.24) shows that all cumulants $\kappa_n$ beyond
$\kappa_2$ vanish; this is a unique property of the Gaussian distribution.
Some other properties are:

$$
\begin{align*}
P(x \text{ in range } \mu \pm \sigma) &= 0.6827, \\
P(x \text{ in range } \mu \pm 0.6745\sigma) &= 0.5, \\
E\left[\left|x - \mu\right|\right] &= \sqrt{2/\pi}\sigma = 0.7979\sigma, \\
\text{half-width at half maximum} &= \sqrt{2 \ln 2} \sigma = 1.177\sigma.
\end{align*}
$$

For a Gaussian with $\mu = 0$ and $\sigma^2 = 1$ (the standard normal)
the cumulative distribution, often written $\Phi(x)$, is related to the
error function erf by

$$
F(x; 0, 1) \equiv \Phi(x) = \frac{1}{2} \left[ 1 + \erf(x/\sqrt{2}) \right] . \quad (39.26)
$$

The error function and standard Gaussian are tabulated in many
references (e.g., Ref. [11, 12]) and are available in software pack-
ages such as ROOT [13]. For a mean $\mu$ and variance $\sigma^2$, replace
$x$ by $(x - \mu)/\sigma$. The probability of $x$ in a given range can be
calculated with Eq. (40.76).

For $x$ and $y$ independent and normally distributed, $z = ax + by$
follows a normal p.d.f. $f(z; a\mu_x + b\mu_y, a^2\sigma_x^2 + b^2\sigma_y^2)$; that is, the
weighted means and variances add.

The Gaussian derives its importance in large part from the cen-
tral limit theorem:

If independent random variables $x_1, \dots, x_n$ are distributed according to any p.d.f. with finite mean and variance, then the sum $y = \sum_{i \neq 1}^n x_i$ will have a p.d.f. that approaches a Gaussian for large $n$. If the p.d.f.s of the $x_i$ are not identical, the theorem still holds under somewhat more restrictive conditions. The mean and variance are given by the sums of corresponding terms from the individual $x_i$. Therefore, the sum of a large number of fluctuations $x_i$ will be distributed as a Gaussian, even if the $x_i$ themselves are not.

For a set of $n$ Gaussian random variables $x$ with means $\mu$ and covariances $V_{ij} = \text{cov}[x_i, x_j]$, the p.d.f. for the one-dimensional Gaussian is generalized to

$$
f(x; \mu, V) = \frac{1}{(2\pi)^{n/2} \sqrt{|V|}} \exp \left[ -\frac{1}{2} (x - \mu)^T V^{-1} (x - \mu) \right], \quad (39.27)
$$

where the determinant $|V|$ must be greater than 0. For diagonal $V$ (independent variables), $f(x; \mu, V)$ is the product of the p.d.f.s of $n$ Gaussian distributions.

For $n = 2$, $f(x; \mu, V)$ is

$$
\begin{equation}
\begin{split}
f(x_1, x_2; \mu_1, \mu_2, \sigma_1, \sigma_2, \rho) = \frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}} \\
\times \exp \left\{ -\frac{1}{2(1-\rho^2)} \left[ \frac{(x_1-\mu_1)^2}{\sigma_1^2} - \frac{2\rho(x_1-\mu_1)(x_2-\mu_2)}{\sigma_1\sigma_2} \right. \right. \\
\left. \left. \qquad + \frac{(x_2-\mu_2)^2}{\sigma_2^2} \right] \right\}.
\end{split}
\tag{39.28}
\end{equation}
$$

The characteristic function for the multivariate Gaussian is

$$
\phi(u; \mu, V) = \exp \left[ i\mu \cdot u - \frac{1}{2} u^T V u \right] . \quad (39.29)
$$

If the components of $\boldsymbol{x}$ are independent, then Eq. (39.29) is the product of the characteristic functions of $n$ Gaussians.

For an $n$-dimensional Gaussian distribution for $\boldsymbol{x}$ with mean $\boldsymbol{\mu}$ and covariance matrix $V$, the marginal distribution for any single $\boldsymbol{x}_i$ is a one-dimensional Gaussian with mean $\boldsymbol{\mu}_i$ and variance $V_{ii}$. The equation $(\boldsymbol{x} - \boldsymbol{a})^T V^{-1} (\boldsymbol{x} - \boldsymbol{a}) = C$, where $C$ is any positive number, defines an $n$-dimensional ellipse centered about $\boldsymbol{a}$. If $\boldsymbol{a}$ is equal to the mean $\boldsymbol{\mu}$, then $C$ is a random variable obeying the $\chi^2$ distribution for $n$ degrees of freedom, which is discussed in the following section. The probability that $\boldsymbol{x}$ lies outside the ellipsoid for a given value of $C$ is given by $1 - F_{\chi^2}(C; n)$, where $F_{\chi^2}$ is the cumulative $\chi^2$ distribution. This may be read from Fig. 40.1. For example, the “s-standard-deviation ellipsoid” occurs at $C = s^2$. For the two-variable case ($n = 2$), the point $\boldsymbol{x}$ lies outside the one-standard-deviation ellipsoid with 61% probability. The use of these ellipsoids as indicators of probable error is described in Sec. 40.4.2.2; the validity of those indicators assumes that $\boldsymbol{\mu}$ and $V$ are correct.

39.4.4 Log-normal distribution

If a random variable $y$ follows a Gaussian distribution with
mean $\mu$ and variance $\sigma^2$, then $x = e^y$ follows a log-normal dis-
tribution, as given in Table 39.1. As a consequence of the cen-
tral limit theorem described in Sec. 39.4.3, the distribution of the
product of a large number of positive random variables approaches
a log-normal. It is bounded below by zero and is thus well suited
for modeling quantities that are intrinsically non-negative such as
an efficiency. One can implement a log-normal model for a ran-
dom variable $x$ by defining $y = \ln x$ so that $y$ follows a Gaussian
distribution.

39.4.5 $\chi^2$ distribution

If $x_1, \dots, x_n$ are independent Gaussian random variables, the
sum $z = \sum_{i=1}^n (x_i - \mu_i)^2 / \sigma_i^2$ follows the $\chi^2$ p.d.f. with $n$ degrees
<table><caption>**Table 39.1:** Some common probability density functions, with corresponding characteristic functions and means and variances. In the Table, $\Gamma(k)$ is the gamma function, equal to $(k-1)!$ when $k$ is an integer; $1F_1$ is the confluent hypergeometric function of the 1st kind [11].</caption><thead><tr><th>Distribution</th><th>Probability density function<br>f (variable; parameters)</th><th>Characteristic<br>function φ(u)</th><th>Mean</th><th>Variance</th></tr></thead><tbody><tr><td>Uniform</td><td>f(x; a, b) = <table><tr><td>1/(b-a)</td><td>a ≤ x ≤ b</td></tr><tr><td>0</td><td>otherwise</td></tr></table></td><td><table><tr><td>e<sup>ibu</sup> - e<sup>iau</sup></td><td>(b-a)iu</td></tr></table></td><td><table><tr><td>a+b</td><td>2</td></tr></table></td><td><table><tr><td>(b-a)<sup>2</sup></td><td>12</td></tr></table></td></tr><tr><td>Binomial</td><td>f(r; N, p) = <table><tr><td>N!</td><td>r!/(N-r)!</td></tr></table> p<sup>r</sup> q<sup>N-r</sup><br>r = 0, 1, 2, ..., N; 0 ≤ p ≤ 1; q = 1 - p</td><td>(q + pe<sup>iu</sup>)<sup>N</sup></td><td>Np</td><td>Npq</td></tr><tr><td>Multinomial</td><td>f(r<sub>1</sub>, ..., r<sub>m</sub>; N, p<sub>1</sub>, ..., p<sub>m</sub>) = <table><tr><td>N!</td><td>r<sub>1</sub>! ... r<sub>m</sub>!</td></tr></table>p<sub>1</sub><sup>r<sub>1</sub></sup> ... p<sub>m</sub><sup>r<sub>m</sub></sup></td><td>(<table><tr><td>Σ<sup>m</sup><sub>k=1</sub> p<sub>k</sub> e<sup>iu<sub>k</sub></sup></td><td>N</td></tr></table>)<sup>N</sup></td><td>E[r<sub>i</sub>] =<br>Np<sub>i</sub></td><td>cov[r<sub>i</sub>, r<sub>j</sub>] =<br>Np<sub>i</sub>(δ<sub>ij</sub> - p<sub>j</sub>)</td></tr><tr><td>Poisson</td><td>f(n; ν) = <table><tr><td>ν<sup>n</sup>e<sup>-ν</sup></td><td>n = 0, 1, 2, ...; ν > 0</td></tr></table></td><td>exp[ν(e<sup>iu</sup> - 1)]</td><td>ν</td><td>ν</td></tr><tr><td>Normal<br>(Gaussian)</td><td>f(x; μ, σ<sup>2</sup>) = <table><tr><td>1</td><td>√2π</td></tr></table> exp(-(x - μ)<sup>2</sup>/2σ<sup>2</sup>)</td><td>exp(iμu - 1/2σ<sup>2</sup>u<sup>2</sup>)</td><td>μ</td><td>σ<sup>2</sup></td></tr><tr><td>Multivariate<br>Gaussian</td><td>f(x; μ, V) = <table><tr><td>1</td><td>(2π)<sup>n/2</sup>√|V|</td></tr></table> × exp[-1/2(x - μ)<sup>T</sup>V<sup>-1</sup>(x - μ)]<br>-∞ < x<sub>j</sub> < ∞; -∞ < μ<sub>j</sub> < ∞; |V| > 0</td><td>exp[iμ · u - 1/2u<sup>T</sup>Vu]</td><td>u</td><td>V<sub>jk</sub></td></tr><tr><td>Log-normal</td><td>f(x; μ, σ<sup>2</sup>) = <table><tr><td>1</td><td>√2π</td></tr></table> 1/2 × exp(-(ln x - μ)<sup>2</sup>/2σ<sup>2</sup>)<br>0 < x < ∞; -∞ < μ < ∞; σ > 0</td><td>—</td><td>exp(μ + σ<sup>2</sup>/2)</td><td>exp(2μ + σ<sup>2</sup>)<br>× [exp(σ<sup>2</sup>) - 1]</td></tr><tr><td>χ<sup>2</sup></td><td>f(z; n) = <table><tr><td>z<sup>n/2-1</sup>e<sup>-z/2</sup></td><td>2n/2Γ(n/2)</td></tr></table>; z ≥ 0</td><td>(1 - 2iu)<sup>-n/2</sup></td><td>n</td><td>2n</td></tr><tr><td>Student's t</td><td>f(t; n) = <table><tr><td>1</td><td>√nπ</td></tr></table> <table><tr><td>Γ[(n+1)/2]</td><td>Γ(n/2)</td></tr></table> (1 + t<sup>2</sup>/n)<sup>-(n+1)/2</sup><br>-∞ < t < ∞; n not required to be integer</td><td>—</td><td>0<br>for n > 1</td><td>n/(n - 2)<br>for n > 2</td></tr><tr><td>Gamma</td><td>f(x; λ, k) = <table><tr><td>x<sup>k-1</sup>λ<sup>k</sup>e<sup>-λx</sup></td><td>Γ(k)</td></tr></table>; 0 ≤ x < ∞; k not required to be integer</td><td>(1 - iu/λ)<sup>-k</sup></td><td>k/λ</td><td>k/λ<sup>2</sup></td></tr><tr><td>Beta</td><td>f(x; α, β) = <table><tr><td>Γ(α+β)</td><td>Γ(α)Γ(β)</td></tr></table> x<sup>α-1</sup>(1 - x)<sup>β-1</sup><br>0 ≤ x ≤ 1</td><td>1F<sub>1</sub>(α; α + β; iu)</td><td>α/α+β</td><td>αβ/(α+β)<sup>2</sup>(α+β+1)</td></tr></tbody></table>

of freedom, which we denote by $\chi^2(n)$. More generally, for $n$ correlated Gaussian variables as components of a vector $X$ with covariance matrix $V$, $z = X^T V^{-1} X$ follows $\chi^2(n)$ as in the previous section. For a set of $z_i$, each of which follows $\chi^2(n_i)$, $\sum z_i$ follows $\chi^2(\sum n_i)$. For large $n$, the $\chi^2$ p.d.f. approaches a Gaussian with a mean and variance given by $\mu = n$ and $\sigma^2 = 2n$, respectively (here the formulae for $\mu$ and $\sigma^2$ are valid for all $n$).

The $\chi^2$ p.d.f. is often used in evaluating the level of compatibility between observed data and a hypothesis for the p.d.f. that the data might follow. This is discussed further in Sec. 40.3.2 on significance tests.

### 39.4.6 Student's t distribution

Suppose that $y$ and $x_1, \dots, x_n$ are independent and Gaussian distributed with mean 0 and variance 1. We then define

$$z = \sum_{i=1}^{n} x_i^2 \quad \text{and} \quad t = \frac{y}{\sqrt{z/n}}. \qquad (39.30)$$

The variable $z$ thus follows a $\chi^2(n)$ distribution. Then $t$ is distributed according to Student's $t$ distribution with $n$ degrees of freedom, $f(t; n)$, given in Table 39.1.

If defined through gamma functions as in Table 39.1, the parameter $n$ is not required to be an integer. As $n \to \infty$, the distribution approaches a Gaussian, and for $n = 1$ it is a Cauchy or Breit-Wigner distribution.

As an example, consider the sample mean $\bar{x} = \sum x_i/n$ and the sample variance $s^2 = \sum (x_i - \bar{x})^2/(n-1)$ for normally distributed $x_i$ with unknown mean $\mu$ and variance $\sigma^2$. The sample mean has a Gaussian distribution with a variance $\sigma^2/n$, so the variable $(\bar{x} - \mu)/\sqrt{\sigma^2/n}$ is normal with mean 0 and variance 1. The quantity $(n-1)s^2/\sigma^2$ is independent of this and follows $\chi^2(n-1)$.

The ratio

$$t = \frac{(\bar{x} - \mu)/\sqrt{\sigma^2/n}}{\sqrt{(n-1)s^2/\sigma^2(n-1)}} = \frac{\bar{x} - \mu}{\sqrt{s^2/n}} \qquad (39.31)$$

is distributed as $f(t; n - 1)$. The unknown variance $\sigma^2$ cancels, and $t$ can be used to test the hypothesis that the true mean is some particular value $\mu$.

### 39.4.7 Gamma distribution

For a process that generates events as a function of $x$ (e.g., space or time) according to a Poisson distribution, the distance in $x$ from an arbitrary starting point (which may be some particular event) to the $k^{th}$ event follows a gamma distribution, $f(x; \lambda, k)$. The Poisson parameter $\lambda$ is $\lambda$ per unit $x$. The special case $k = 1$ (i.e., $f(x; \lambda, 1) = \lambda e^{-\lambda x}$) is called the exponential distribution. A sum of $k'$ exponential random variables $x_i$ is distributed as $f(\sum x_i; \lambda, k')$.

The parameter $k$ is not required to be an integer. For $\lambda = 1/2$ and $k = n/2$, the gamma distribution reduces to the $\chi^2(n)$ distribution.

### 39.4.8 Beta distribution

The beta distribution describes a continuous random variable $x$ in the interval $[0, 1]$. By scaling and translation one can easily generalize it to have arbitrary endpoints. In Bayesian inference about the parameter $p$ of a binomial process, if the prior p.d.f. is a beta distribution $f(p; \alpha, \beta)$ then the observation of $r$ successes out of $N$ trials gives a posterior beta distribution $f(p; r + \alpha, N - r + \beta)$ (Bayesian methods are discussed further in Sec. 40). The uniform distribution is a beta distribution with $\alpha = \beta = 1$.

### References

[1] H. Cramér, *Mathematical Methods of Statistics*, (Princeton Univ. Press, New Jersey, 1958).
[2] A. Stuart and J.K. Ord, *Kendall's Advanced Theory of Statistics*, Vol. 1 *Distribution Theory* 6th Ed., (Halsted Press, New York, 1994), and earlier editions by Kendall and Stuart.

[3] F.E. James, *Statistical Methods in Experimental Physics*, 2nd Ed., (World Scientific, Singapore, 2006).

[4] L. Lyons, *Statistics for Nuclear and Particle Physicists*, (Cambridge University Press, New York, 1986).

[5] B.R. Roe, *Probability and Statistics in Experimental Physics*, 2nd Ed., (Springer, New York, 2001).

[6] R.J. Barlow, *Statistics: A Guide to the Use of Statistical Methods in the Physical Sciences*, (John Wiley, New York, 1989).

[7] S. Brandt, *Data Analysis*, 3rd Ed., (Springer, New York, 1999).

[8] G. Cowan, *Statistical Data Analysis*, (Oxford University Press, Oxford, 1998).

[9] A.N. Kolmogorov, *Grundbegriffe der Wahrscheinlichkeitsrechnung*, (Springer, Berlin, 1933); *Foundations of the Theory of Probability*, 2nd Ed., (Chelsea, New York 1956).

[10] Ch. Walck, *Hand-book on Statistical Distributions for Experimentalists*, University of Stockholm Internal Report SUF-PFY/96-01, available from www.physto.se-walck.

[11] M. Abramowitz and I. Stegun, eds., *Handbook of Mathematical Functions*, (Dover, New York, 1972).

[12] F.W.J. Olver et al., eds., NIST Handbook of Mathematical Functions, (Cambridge University Press, 2010); A companion Digital Library of Mathematical Functions is available at dlmf.nist.gov.

[13] R. Brun and F. Rademakers, Nucl. Instrum. Meth. **A389**, 81 (1997); See also root.cern.ch.
# 40. Statistics

Revised August 2023 by G. Cowan (RHUL).

This chapter gives an overview of statistical methods used in high-energy physics. In statistics, we are interested in using a given sample of data to make inferences about a probabilistic model, e.g., to assess the model's validity or to determine the values of its parameters. There are two main approaches to statistical inference, which we may call frequentist and Bayesian.

In frequentist statistics, probability is interpreted as the limiting frequency of the outcome of a repeatable experiment. The most important tools in this framework are parameter estimation, covered in Section 40.2, statistical tests, discussed in Section 40.3, and confidence intervals, which are constructed so as to cover the true value of a parameter with a specified probability, as described in Section 40.4.2. Note that in frequentist statistics one does not define a probability for a hypothesis or for the value of a parameter.

In Bayesian statistics, the subjective interpretation of probability is used to quantify one's degree of belief in a hypothesis. This allows one to define a probability density function (p.d.f.) for a parameter, which reflects one's knowledge about where its true value lies.

Bayesian methods provide a natural means to include additional information, which in general may be subjective; in fact they require prior probabilities for the hypotheses (or parameters) in question, i.e., the degree of belief about the parameters' values, before carrying out the measurement. Using Bayes' theorem (Eq. (39.4)), the prior degree of belief is updated by the data from the experiment. Bayesian methods for interval estimation are discussed in Sections 40.4.1 and 40.4.2.4.

For many inference problems, the frequentist and Bayesian approaches give similar numerical values, even though they answer different questions and are based on fundamentally different interpretations of probability. In some important cases, however, the two approaches may yield very different results. For a discussion of Bayesian vs. non-Bayesian methods, see references written by a statistician [1], by a physicist [2], or the detailed comparison in Ref. [3].

## 40.1 Fundamental concepts

Consider an experiment whose outcome is characterized by one or more data values, which we can write as a vector $x$. A hypothesis $H$ is a statement about the probability for the data, often written $P(x|H)$. (We will usually use a capital letter for a probability and lower case for a probability density. Often the term p.d.f. is used loosely to refer to either a probability or a probability density.) This could, for example, define completely the p.d.f. for the data (a simple hypothesis), or it could specify only the functional form of the p.d.f., with the values of one or more parameters not determined (a composite hypothesis).

If the probability $P(x|H)$ for data $x$ is regarded as a function of the hypothesis $H$, then it is called the likelihood of $H$, usually written $L(H)$. Often the hypothesis is characterized by one or more parameters $\theta$, in which case $L(\theta) = P(x|\theta)$ is called the likelihood function.

In some cases one can obtain at least approximate frequentist results using the likelihood evaluated only with the data obtained, for example, when determining confidence regions with Eq. (40.79). In general, however, the frequentist approach requires a full specification of the probability model $P(x|H)$ both as a function of the data $x$ and hypothesis $H$.

In the Bayesian approach, inference is based on the posterior probability for $H$ given the data $x$, which represents one's degree of belief that $H$ is true given the data. This is obtained from Bayes' theorem (39.4), which can be written

$$P(H|x) = \frac{P(x|H)\pi(H)}{\int P(x|H')\pi(H') dH'} \quad (40.1)$$

Here $P(x|H)$ is the likelihood for $H$, which depends only on the data actually obtained. The quantity $\pi(H)$ is the prior probability for $H$, which represents one's degree of belief for $H$ before carrying out the measurement. The integral in the denominator (or sum, for discrete hypotheses) serves as a normalization factor. If $H$

is characterized by a continuous parameter $\theta$ then the posterior probability is a p.d.f. $p(\theta|x)$. Note that the likelihood function itself is not a p.d.f. for $\theta$.

## 40.2 Parameter estimation

Here we review point estimation of parameters, first with an overview of the frequentist approach and its two most important methods, maximum likelihood and least squares, treated in Sections 40.2.2 and 40.2.3. The Bayesian approach is outlined in Sec. 40.2.6.

An estimator $\hat{\theta}$ (written with a hat) is a function of the data used to estimate the value of the parameter $\theta$. Sometimes the word 'estimate' is used to denote the value of the estimator when evaluated with given data. There is no fundamental rule dictating how an estimator must be constructed. One tries, therefore, to choose that estimator which has the best properties. The most important of these are (a) consistency, (b) bias, (c) efficiency, and (d) robustness.

(a) An estimator is said to be *consistent* if the estimate $\hat{\theta}$ converges in probability (see Ref. [3]) to the true value $\theta$ as the amount of data increases. This property is so important that it is possessed by all commonly used estimators.

(b) The bias, $b = E[\hat{\theta}] - \theta$, is the difference between the expectation value of the estimator and the true value of the parameter. The expectation value is taken over a hypothetical set of similar experiments in which $\hat{\theta}$ is constructed in the same way. When $b=0$, the estimator is said to be unbiased. The bias depends on the chosen metric, i.e., if $\hat{\theta}$ is an unbiased estimator of $\theta$, then $\hat{\theta}^2$ is not in general an unbiased estimator for $\theta^2$.

(c) Efficiency is the ratio of the minimum possible variance for any estimator of $\theta$ to the variance $V[\hat{\theta}]$ of the estimator $\hat{\theta}$. For the case of a single parameter, under rather general conditions the minimum variance is given by the Rao-Cramér-Fréchet bound,

$$\sigma_{\min}^2 = \left(1 + \frac{\partial b}{\partial \theta}\right)^2 / I(\theta), \quad (40.2)$$

where

$$I(\theta) = E \left[ \left( \frac{\partial \ln L}{\partial \theta} \right)^2 \right] = -E \left[ \frac{\partial^2 \ln L}{\partial \theta^2} \right] \quad (40.3)$$

is the Fisher information, $L$ is the likelihood, and the operator $E[]$ in (40.3) is the expectation value with respect to the data. For the final equality to hold, the range of allowed data values must not depend on $\theta$.

The mean-squared error,

$$\text{MSE} = E[(\hat{\theta} - \theta)^2] = V[\hat{\theta}] + b^2, \quad (40.4)$$

is a measure of an estimator's quality which combines bias and variance.

(d) Robustness is the property of being insensitive to departures from assumptions about the p.d.f., e.g., owing to uncertainties in the distribution's tails.

It is not in general possible to optimize simultaneously for all the measures of estimator quality described above. For some common estimators, the properties above are known exactly. More generally, it is possible to evaluate them by Monte Carlo simulation. Note that they will in general depend on the unknown $\theta$.

### 40.2.1 Estimators for mean, variance, and median

Suppose we have a set of $n$ independent measurements, $x_1, \dots, x_n$, each assumed to follow a p.d.f. with unknown mean $\mu$ and unknown variance $\sigma^2$ (the measurements do not necessarily have to follow a Gaussian distribution). Then

$$\hat{\mu} = \frac{1}{n} \sum_{i=1}^{n} x_i \quad (40.5)$$

$$\hat{\sigma}^2 = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \hat{\mu})^2 \quad (40.6)$$
are unbiased estimators of $\mu$ and $\sigma^2$. The variance of $\hat{\mu}$ is $\sigma^2/n$ and the variance of $\hat{\sigma}^2$ is

$$V[\hat{\sigma}^2] = \frac{1}{n} \left( m_4 - \frac{n-3}{n-1} \sigma^4 \right), \quad (40.7)$$

where $m_4$ is the 4th central moment of $x$ (see Eq. (39.8)). For Gaussian distributed $x_i$, this becomes $2\sigma^4/(n-1)$ for any $n \ge 2$, and for large $n$ the standard deviation of $\hat{\sigma}$ is $\sigma/\sqrt{2n}$. For any $n$ and Gaussian $x_i$, $\hat{\mu}$ is an efficient estimator for $\mu$, and the estimators $\hat{\mu}$ and $\hat{\sigma}^2$ are uncorrelated. Otherwise the arithmetic mean (40.5) is not necessarily the most efficient estimator; this is discussed further in Sec. 8.7 of Ref. [4].

If $\sigma^2$ is known, it does not improve the estimate $\hat{\mu}$, as can be seen from Eq. (40.5); however, if $\mu$ is known, one can substitute it for $\hat{\mu}$ in Eq. (40.6) and replace $n-1$ by $n$ to obtain an estimator of $\sigma^2$ still with zero bias but smaller variance. If the $x_i$ have different, known variances $\sigma_i^2$, then the weighted average

$$\hat{\mu} = \frac{1}{w} \sum_{i=1}^{n} w_i x_i, \quad (40.8)$$

where $w_i = 1/\sigma_i^2$ and $w = \sum_i w_i$, is an unbiased estimator for $\mu$ with a smaller variance than an unweighted average. The standard deviation of $\hat{\mu}$ is $1/\sqrt{w}$.

As an estimator for the median $x_{\text{med}}$, one can use the value $\hat{x}_{\text{med}}$ such that half the $x_i$ are below and half above (the sample median). If there are an even number of observations and the sample median lies between two observed values, the estimator is set by convention to their arithmetic average. If the p.d.f. of $x$ has the form $f(x - \mu)$ and $\mu$ is both mean and median, then for large $n$ the variance of the sample median approaches $1/[4nf^2(0)]$, provided $f(0) > 0$ [5]. Although estimating the median can often be more difficult computationally than the mean, the resulting estimator is generally more robust, as it is insensitive to the exact shape of the tails of a distribution.

### 40.2.2 The method of maximum likelihood

Suppose we have a set of measured quantities $x$ and the likelihood $L(\theta) = P(x|\theta)$ for a set of parameters $\theta = (\theta_1, \dots, \theta_M)$. The maximum likelihood estimators (MLEs) for $\theta$ are defined as the values that give the maximum of $L$. Because of the properties of the logarithm, it is usually easier to work with $\ln L$, and since both are maximized for the same parameter values $\theta$, the MLEs can be found by solving the likelihood equations,

$$\frac{\partial \ln L}{\partial \theta_i} = 0, \quad i = 1, \dots, M. \quad (40.9)$$

Often the solution must be found numerically. Maximum likelihood estimators are important because they are asymptotically (i.e., for large data samples) unbiased, efficient and have a Gaussian sampling distribution under quite general conditions, and the method has a wide range of applicability.

In general the likelihood function is obtained from the probability of the data under assumption of the parameters. An important special case is when the data consist of i.i.d. (independent and identically distributed) values. Here one has a set of $n$ statistically independent quantities $x = (x_1, \dots, x_n)$, where each component follows the same p.d.f. $f(x; \theta)$. In this case the joint p.d.f. of the data sample factorizes and the likelihood function is

$$L(\theta) = \prod_{i=1}^{n} f(x_i; \theta). \quad (40.10)$$

Here the number $n$ of observations (usually individual “events” in particle physics) is regarded as fixed. If, however, the probability to observe $n$ events itself depends on the parameters $\theta$, then this dependence should be included in the likelihood. For example, if $n$ follows a Poisson distribution with mean $\mu$ and the independent $x$ values all follow $f(x; \theta)$, then the likelihood becomes

$$L(\theta) = \frac{\mu^n}{n!} e^{-\mu} \prod_{i=1}^{n} f(x_i; \theta). \quad (40.11)$$

Equation (40.11) is often called the extended likelihood (see, e.g., Refs. [6–8]). If $\mu$ is given as a function of $\theta$, then including the probability for $n$ given $\theta$ in the likelihood provides additional information about the parameters. This therefore leads to a reduction in their statistical uncertainties and in general changes their estimated values.

In evaluating the likelihood function, it is important that any normalization factors in the p.d.f. that involve $\theta$ be included. However, we will only be interested in the maximum of $L$ and in ratios of $L$ at different values of the parameters; hence any multiplicative factors that do not involve the parameters that we want to estimate may be dropped, including factors that depend on the data but not on $\theta$.

Under a one-to-one change of parameters from $\theta$ to $\eta$, the MLEs $\hat{\theta}$ transform to $\eta(\hat{\theta})$. That is, the maximum-likelihood solution is invariant under change of parameter. However, other properties of MLEs, in particular the bias, are not invariant under change of parameter.

The inverse $V^{-1}$ of the covariance matrix $V_{ij} = \text{cov}[\hat{\theta}_i, \hat{\theta}_j]$ for a set of MLEs can be estimated by using

$$(\hat{V}^{-1})_{ij} = - \frac{\partial^2 \ln L}{\partial \theta_i \partial \theta_j} \Big|_{\hat{\theta}}. \quad (40.12)$$

Equation (40.12) holds for a sufficiently large data sample (or for a small sample only in special cases, e.g., where the means of Gaussian distributed data are linear functions of the parameters); otherwise it can result in a misestimation of the variances. Under the conditions where the equation is valid, $L$ has a Gaussian form and $\ln L$ is (hyper)parabolic. In this case, $s$ times the standard deviations $\sigma_i$ of the estimators for the parameters can be obtained from the hypersurface defined by the $\theta$ such that

$$\ln L(\theta) = \ln L_{\text{max}} - s^2/2, \quad (40.13)$$

where $\ln L_{\text{max}}$ is the value of $\ln L$ at the solution point (compare with Eq. (40.79)). The minimum and maximum values of $\theta_i$ on the hypersurface then give an approximate $s$-standard deviation confidence interval for $\theta_i$ (see Section 40.4.2.2).

#### 40.2.2.1 Maximum likelihood with binned data

If the total number of data values $x_i$, ($i = 1, \dots, n_{\text{tot}}$), is small, the unbinned maximum-likelihood method, i.e., use of Equation (40.10) (or (40.11) for extended maximum likelihood), is preferred since binning can only result in a loss of information, and hence larger statistical errors for the parameter estimates. If the sample is large, it can be convenient to bin the values in a histogram with $N$ bins, so that one obtains a vector of data $n = (n_1, \dots, n_N)$ with expectation values $\mu = E[n]$ and probabilities $f(n; \mu)$. Suppose the mean values $\mu$ can be determined as a function of a set of parameters $\theta$. Then one may maximize the likelihood function based on the contents of the bins.

As mentioned in Sec. 40.2.2, the total number of events $n_{\text{tot}} = \sum_i n_i$ can be regarded either as fixed or as a random variable. If it is fixed, the histogram follows a multinomial distribution,

$$f_{\text{M}}(n; \theta) = \frac{n_{\text{tot}}!}{n_1! \cdots n_N!} p_1^{n_1} \cdots p_N^{n_N}, \quad (40.14)$$

where we assume the probabilities $p_i$ are given functions of the parameters $\theta$. The distribution can be written equivalently in terms of the expected number of events in each bin, $\mu_i = n_{\text{tot}}p_i$. If the $n_i$ are regarded as independent and Poisson distributed, then the data are instead described by a product of Poisson probabilities,

$$f_{\text{P}}(n; \theta) = \prod_{i=1}^{N} \frac{\mu_i^{n_i}}{n_i!} e^{-\mu_i}, \quad (40.15)$$

where the mean values $\mu_i$ are given functions of $\theta$. The total number of events $n_{\text{tot}}$ thus follows a Poisson distribution with mean $\mu_{\text{tot}} = \sum_i \mu_i$.

When using maximum likelihood with binned data, one can find the estimators and at the same time obtain a statistic usable for a test of goodness of fit (see Sec. 40.3.3.1). For independent Poisson distributed data, maximizing the likelihood
$L(\theta) = f_P(n; \theta)$ is equivalent to maximizing the likelihood ratio $\lambda(\theta) = f_P(n; \theta)/f_{sat}(n; \hat{\mu})$, where in the denominator $f_{sat}(n; \mu)$ is the corresponding model with an adjustable parameter for each bin, $\mu = (\mu_1, \dots, \mu_N)$, with estimators $\hat{\mu} = (n_1, \dots, n_N)$ (called the “saturated model”). Equivalently one often minimizes the quantity $-2 \ln \lambda(\theta)$. For independent Poisson distributed $n_i$ this is [9]

$$-2 \ln \lambda(\theta) = 2 \sum_{i=1}^{N} \left[ \mu_i(\theta) - n_i + n_i \ln \frac{n_i}{\mu_i(\theta)} \right], \quad (40.16)$$

where for bins with $n_i = 0$, the last term in (40.16) is zero. The expression (40.16) without the terms $\mu_i - n_i$ also gives $-2 \ln \lambda(\theta)$ for multinomially distributed $n_i$, i.e., when the total number of entries is regarded as fixed. In the limit of zero bin width, minimizing (40.16) is equivalent to maximizing the unbinned extended likelihood function (40.11); in the corresponding multinomial case without the $\mu_i - n_i$ terms one obtains Eq. (40.10).

A smaller value of $-2 \ln \lambda(\hat{\theta})$ corresponds to better agreement between the data and the hypothesized form of $\mu(\theta)$. The value of $-2 \ln \lambda(\hat{\theta})$ can thus be translated into a p-value as a measure of goodness-of-fit, as described in Sec. 40.3.3.1. Assuming the model is correct, then according to Wilks' theorem [10], for sufficiently large $\mu_i$ and provided certain regularity conditions are met, the minimum of $-2 \ln \lambda$ as defined by Eq. (40.16) follows a $\chi^2$ distribution (see, e.g., Ref. [9]). If there are $N$ bins and $M$ fitted parameters, then the number of degrees of freedom for the $\chi^2$ distribution is $N - M$ if the data are treated as Poisson-distributed, and $N - M - 1$ if the $n_i$ are multinomially distributed.

Suppose the $n_i$ are Poisson-distributed and the overall normalization $\mu_{tot} = \sum_i \mu_i$ is taken as an adjustable parameter, so that $\mu_i = \mu_{tot}p_i(\theta)$, where the probability to be in the $i$th bin, $p_i(\theta)$, does not depend on $\mu_{tot}$. Then by minimizing Eq. (40.16), one obtains that the area under the fitted function is equal to the sum of the histogram contents, i.e., $\sum_i \hat{\mu}_i = \sum_i n_i$. This is a property not possessed by the estimators from the method of least squares (see, e.g., Sec. 40.2.3 and Ref. [8]).

### 40.2.2.2 Frequentist treatment of nuisance parameters

Suppose we want to determine the values of parameters $\theta$ using a set of measurements $x$ described by a probability model $P(x|\theta)$. In general the model is not perfect, which is to say it cannot provide an accurate description of the data even at the most optimal point of its parameter space. As a result, the estimated parameters can have a systematic bias.

One can improve the model by including in it additional parameters. That is, $P(x|\theta)$ is replaced by a more general model $P(x|\theta, \nu)$, which depends on parameters of interest $\theta$ and nuisance parameters $\nu$. The additional parameters are not of intrinsic interest but must be included for the model to be sufficiently accurate for some point in the enlarged parameter space.

Although including additional parameters may eliminate or at least reduce the effect of systematic uncertainties, their presence will result in increased statistical uncertainties for the parameters of interest. This occurs because the estimators for the nuisance parameters and those of interest will in general be correlated, which results in an enlargement of the contour defined by Eq. (40.13).

To reduce the impact of the nuisance parameters one often tries to constrain their values by means of control or calibration measurements, say, having data $y$. For example, some components of $y$ could represent estimates of the nuisance parameters, often from separate experiments. Suppose the measurements $y$ are statistically independent from $x$ and are described by a model $P(y|\nu)$. The joint model for both $x$ and $y$ is in this case therefore the product of the probabilities for $x$ and $y$, and thus the likelihood function for the full set of parameters is

$$L(\theta, \nu) = P(x|\theta, \nu)P(y|\nu). \quad (40.17)$$

Note that in this case if one wants to simulate the experiment by means of Monte Carlo, both the primary and control mea-

surements, $x$ and $y$, must be generated for each repetition under assumption of fixed values for the parameters $\theta$ and $\nu$.

Using all of the parameters $(\theta, \nu)$ in Eq. (40.13) to find the statistical errors in the parameters of interest $\theta$ is equivalent to using the profile likelihood, which depends only on $\theta$. It is defined as

$$L_p(\theta) = L(\theta, \hat{\nu}(\theta)), \quad (40.18)$$

where the double-hat notation indicates the profiled values of the parameters $\nu$, defined as the values that maximize $L$ for the specified $\theta$. The profile likelihood is discussed further in Section 40.3.2.1 in connection with hypothesis tests.

### 40.2.3 The method of least squares

The method of least squares (LS) coincides with the method of maximum likelihood in the following special case. Consider a set of $N$ independent measurements $y_i$ at known points $x_i$. The measurement $y_i$ is assumed to be Gaussian distributed with mean $\mu(x_i; \theta)$ and known variance $\sigma_i^2$. The goal is to construct estimators for the unknown parameters $\theta$. The log-likelihood function contains the sum of squares

$$\chi^2(\theta) = -2 \ln L(\theta) + \text{constant} = \sum_{i=1}^{N} \frac{(y_i - \mu(x_i; \theta))^2}{\sigma_i^2} . \quad (40.19)$$

The parameter values that maximize $L$ are the same as those which minimize $\chi^2$.

The minimum of the chi-square function in Equation (40.19) defines the least-squares estimators $\hat{\theta}$ for the more general case where $y_i$ are not Gaussian distributed as long as they are independent. If they are not independent but rather have a covariance matrix $V_{ij} = \text{cov}[y_i, y_j]$, then the LS estimators are determined by the minimum of

$$\chi^2(\theta) = (y - \mu(\theta))^T V^{-1} (y - \mu(\theta)), \quad (40.20)$$

where $y = (y_1, \dots, y_N)$ is the (column) vector of measurements, $\mu(\theta)$ is the corresponding vector of predicted values, and the superscript $T$ denotes the transpose. If the $y_i$ are not Gaussian distributed, then the least-squares and maximum-likelihood estimators will not in general coincide.

Often one further restricts the problem to the case where $\mu(x_i; \theta)$ is a linear function of the parameters, i.e.,

$$\mu(x_i; \theta) = \sum_{j=1}^{m} \theta_j h_j(x_i). \quad (40.21)$$

Here the $h_j(x)$ are $m$ linearly independent functions, e.g., $1, x, x^2, \dots, x^{m-1}$ or Legendre polynomials. We require $m < N$ and at least $m$ of the $x_i$ must be distinct.

Minimizing $\chi^2$ in this case with $m$ parameters reduces to solving a system of $m$ linear equations. Defining $H_{ij} = h_j(x_i)$ and minimizing $\chi^2$ by setting its derivatives with respect to the $\theta_i$ equal to zero gives the LS estimators,

$$\hat{\theta} = (H^T V^{-1} H)^{-1} H^T V^{-1} y \equiv Dy. \quad (40.22)$$

The covariance matrix for the estimators $U_{ij} = \text{cov}[\hat{\theta}_i, \hat{\theta}_j]$ is given by

$$U = DVD^T = (H^T V^{-1} H)^{-1}, \quad (40.23)$$

or equivalently, its inverse $U^{-1}$ can be found from

$$(U^{-1})_{ij} = \frac{1}{2} \frac{\partial^2 \chi^2}{\partial \theta_i \partial \theta_j} \Big|_{\theta = \hat{\theta}} = \sum_{k,l=1}^{N} h_i(x_k) (V^{-1})_{kl} h_j(x_l). \quad (40.24)$$

The LS estimators can also be found from the expression

$$\hat{\theta} = U g, \quad (40.25)$$

where the vector $g$ is defined by

$$g_i = \sum_{j,k=1}^{N} y_j h_i(x_k) (V^{-1})_{jk}. \quad (40.26)$$
For the case of uncorrelated $y_i$, for example, one can use (40.25) with

$$ (U^{-1})_{ij} = \sum_{k=1}^{N} \frac{h_i(x_k)h_j(x_k)}{\sigma_k^2}, \quad (40.27) $$

$$ g_i = \sum_{k=1}^{N} \frac{y_k h_i(x_k)}{\sigma_k^2}. \quad (40.28) $$

Expanding $\chi^2(\theta)$ about $\hat{\theta}$, one finds that the contour in parameter space defined by

$$ \chi^2(\theta) = \chi^2(\hat{\theta}) + 1 = \chi^2_{\min} + 1 \quad (40.29) $$

has tangent planes located at plus-or-minus-one standard deviation $\sigma_{\hat{\theta}}$ from the LS estimates $\hat{\theta}$ (the relation is approximate if the fit function $\mu(x; \theta)$ is nonlinear in the parameters).

In constructing the quantity $\chi^2(\theta)$ one requires the variances or, in the case of correlated measurements, the covariance matrix. Often these quantities are not known a priori and must be estimated from the data. In this case the least-squares and maximum-likelihood methods are no longer exactly equivalent even for Gaussian distributed measurements. An important example is where the measured value $y_i$ represents the event count in a histogram bin. If, for example, $y_i$ represents a Poisson variable, for which the variance is equal to the mean, then one can either estimate the variance from the predicted value, $\mu(x_i; \theta)$, or from the observed number itself, $y_i$. In the first option, the variances become functions of the parameters, and as a result the estimators may need to be found numerically. The second option can be undefined if $y_i$ is zero, and for small $y_i$, the variance will be poorly estimated. In either case, one should constrain the normalization of the fitted curve to the correct value, i.e., one should determine the area under the fitted curve directly from the number of entries in the histogram (see Ref. [8], Section 7.4). As noted in Sec. 40.2.2.1, this issue is avoided when using the method of extended maximum likelihood with binned data by minimizing Eq. (40.16). In that case if the expected number of events $\mu_{tot}$ does not depend on the other fitted parameters $\theta$, then its extended MLE is equal to the observed total number of events.

As the minimum value of the $\chi^2$ represents the level of agreement between the measurements and the fitted function, it can be used for assessing the goodness-of-fit; this is discussed further in Section 40.3.3.1.

### 40.2.4 Parameter estimation with constraints

In some applications one is interested in using a set of measured quantities $y = (y_1, \dots, y_N)$ to estimate parameters $\theta = (\theta_1, \dots, \theta_M)$ subject to a number of constraints. For example, one may have measured coordinates from two tracks, and the goal is to estimate their momentum vectors subject to the constraint that the tracks have a common vertex. The parameters can also include momenta of undetected particles such as neutrons, as long as the constraints from conservation of energy and momentum and from known masses of particles involved in the reaction chain provide enough information for these quantities to be inferred.

A set of $K$ constraints can be given in the form of equations

$$ c_k(\theta) = 0, \quad k = 1, \dots, K. \quad (40.30) $$

In some problems it may be possible to define a new set of parameters $\eta = (\eta_1, \dots, \eta_L)$ with $L = M - K$ such that every point in $\eta$-space automatically satisfies the constraints. If this is possible then the problem reduces to one of estimating $\eta$ with, e.g., maximum likelihood or least squares and then transforming the estimators back into $\theta$-space.

In many cases it may be difficult or impossible to find an appropriate transformation $\eta(\theta)$. Suppose that the parameters are determined by minimizing an objective function such as $\chi^2(\theta)$ in the method of least squares. Here one may enforce the constraints

by finding the stationary points of the Lagrange function

$$ \mathcal{L}(\theta, \lambda, y) = \chi^2(\theta, y) + \sum_{k=1}^{K} \lambda_k c_k(\theta) \quad (40.31) $$

with respect to both the parameters $\theta$ and a set of Lagrange multipliers $\lambda = (\lambda_1, \dots, \lambda_K)$. Combining the parameters and Lagrange multipliers into an $(M + K)$-component vector $\gamma = (\theta_1, \dots, \theta_M, \lambda_1, \dots, \lambda_K)$, the solutions for $\gamma$, i.e., the estimators $\hat{\gamma}$, are found (e.g., numerically) from the system of equations

$$ F_i(\gamma, y) \equiv \frac{\partial \mathcal{L}}{\partial \gamma_i} = 0, \quad i = 1, \dots, M + K. \quad (40.32) $$

To obtain the covariance matrix of the estimated parameters one can find solutions $\tilde{\gamma}$ corresponding to the expectation values of the data $\langle y \rangle$ and expand $F_i(\tilde{\gamma}, y)$ to first order about these values. This gives (see, e.g., Sec. 11.6 of Ref. [8]) linearized approximations for the estimators, $\hat{\gamma}(y) \approx \tilde{\gamma} + C(y - \langle y \rangle)$, where the matrix $C = -A^{-1}B$, and $A$ and $B$ are given by

$$ A_{ij} = \left[ \frac{\partial F_i}{\partial \gamma_j} \right]_{\tilde{\gamma}, \langle y \rangle} \quad \text{and} \quad B_{ij} = \left[ \frac{\partial F_i}{\partial y_j} \right]_{\tilde{\gamma}, \langle y \rangle} \quad (40.33) $$

In practice the values $\langle y \rangle$ and corresponding solutions $\tilde{\gamma}$ are estimated using the data from the actual measurement. Using this approximation for $\hat{\gamma}(y)$, one can find the covariance matrix $U_{ij} = \text{cov}[\hat{\gamma}_i, \hat{\gamma}_j]$ of the estimators for the $\gamma_i$ in terms of that of the data $V_{ij} = \text{cov}[y_i, y_j]$ using error propagation (cf. Eqs. (39.17) and (39.18)),

$$ U = CVC^T. \quad (40.34) $$

The upper-left $M \times M$ block of the matrix $U$ gives the covariance matrix for the estimated parameters $\text{cov}[\hat{\theta}_i, \hat{\theta}_j]$. One can show for linear constraints that $\text{cov}[\hat{\theta}_i, \hat{\theta}_j]$ is also given by the upper-left $M \times M$ block of $2A^{-1}$. If the parameters are estimated using the method of least squares, then the number of degrees of freedom for the distribution of the minimized $\chi^2$ is increased by the number of constraints, i.e., it becomes $N - M + K$. Further details can be found in, e.g., Ch. 8 of Ref. [4] and Ch. 7 of Ref. [11].

### 40.2.5 Unfolding

An important class of parameter estimation problem involves measurement of the differential distribution of a kinematic variable in the form of a histogram with $N$ bins. The data thus consist of the vector of measured data values $n = (n_1, \dots, n_N)$, with expectation values $\nu_i = E[n_i]$. The $n_i$ are usually independent and often modeled as Poisson distributed, from which it follows that the maximum-likelihood estimators are $\hat{\nu}_i = n_i$ for all $i$.

Because of the limited acceptance and resolution of the experiment, however, the measured values of the kinematic variable in question differ in general from their true values, and as a consequence the form of the data histogram is distorted relative to what would be obtained with perfect resolution. The desired parameters are not, therefore, the $\nu_i$ but rather one wants to estimate the expected number of entries in a given bin that would be found with a perfect detector. We call this the “true histogram” and denote it with $\mu = (\mu_1, \dots, \mu_M)$, where the number of bins $M$ is not required to equal that of the observed histogram. The $\mu_i$ are related to the expected numbers of events in the observed histogram by

$$ \nu_i = \sum_{j=1}^{M} R_{ij} \mu_j + \beta_i, \quad (40.35) $$

where $\beta_i$ here represents the expected number of events in bin $i$ due to background processes and the $N \times M$ response matrix $R_{ij}$ gives the probability for an event to be observed in bin $i$ of $\nu$ given that the true value of the variable was in bin $j$ of $\mu$. For purposes of this discussion let us suppose that the response matrix and expected background values are known. There are two main approaches to this type of problem, which we can call unfolding and folding.
In unfolding, one treats the true histogram $\mu$ as the parameters of interest. The result is thus given by estimators $\hat{\mu}$ and the corresponding covariance matrix $U_{ij} = \text{cov}[\hat{\mu}_i, \hat{\mu}_j]$. Provided the response matrix can be inverted, the maximum-likelihood solution is easily found as $\hat{\mu} = R^{-1}(n - \beta)$. If the response matrix allows for significant migration of events between bins, then the variances of these estimators can be very large, sometimes to the point where the $\hat{\mu}$ bear essentially no resemblance to the true $\mu$. In such cases the estimators can be found by maximizing a linear combination of the log-likelihood and a regularization function that imposes some degree of smoothness on the unfolded distribution. In achieving a reduction in variance one inevitably introduces some bias into the estimators.

In the approach of folding, by contrast, to test a given model prediction for the true histogram $\mu$ it is first “folded” with the response matrix and corrected for expected background to give the corresponding $\nu$ according to Eq. (40.35), and these are compared to the corresponding $n$ using, e.g., a likelihood ratio as shown in Sec. 40.2.2.1. In this way, one avoids any bias due to regularization.

To account for systematic uncertainties, the response matrix and background values may not be expressed as constants but rather as functions of nuisance parameters $\theta$. These may be constrained by auxiliary measurements $u$ with p.d.f. $p(u|\theta)$ so that the full likelihood becomes (as in Eq. (40.17) but with different notation)

$$L(\mu, \theta) = P(n|\nu(\mu, \theta))p(u|\theta). \quad (40.36)$$

Here the mean values $\nu = R(\theta)\mu + \beta(\theta)$ now depend on both the true histogram parameters $\mu$ as well as the nuisance parameters $\theta$. Thus to record enough information for a future analysis using folding one must provide all of the ingredients used above: the primary data $n$, the auxiliary measurements $u$, and also the response matrix $R(\theta)$ and expected backgrounds $\beta(\theta)$ as functions of the nuisance parameters. If unfolding is used, then future model tests or comparisons with other experiments can be carried out directly using the estimators $\hat{\mu}$ and their covariance matrix.

If several distributions are unfolded, then to combine these in a test of a given model one should know how estimators for bins of different distributions are correlated, as can arise, e.g., through common systematic effects. If only the unfolded distributions and their separate covariance matrices are reported, however, then information on such correlations is not retained. In folding, one can include information on correlated systematic effects if the ingredients ($R$ and $\beta$) are known in terms of nuisance parameters that are common to different distributions.

In unfolding, the estimators for $\mu$ can be constructed either by maximizing the log-likelihood using Eq. (40.36) or in the case where regularization is required one can maximize $\varphi(\mu, \theta) = \ln L(\mu, \theta) + \tau S(\mu)$, where the regularization function $S(\mu)$ reflects the smoothness of the true histogram and serves to reduce the variances of the estimators. Possible functions are based on the mean squared second derivative (a commonly used type of Tikhonov regularization) or the entropy of the true histogram. The parameter $\tau$ fixes the relative weighting of the log-likelihood and the regularization function and thus determines the balance between the bias and variance of the resulting estimators. Further discussion of the unfolding problem including methods for choosing the regularization function and parameter, as well as techniques that employ other types of regularization such as the iterative Bayes (Richardson-Lucy) method, can be found in Refs. [8, 11–13] and references therein.

### 40.2.6 The Bayesian approach

In the frequentist methods discussed above, probability is associated only with data, not with the value of a parameter. This is no longer the case in Bayesian statistics, however, which we introduce in this section. For general introductions to Bayesian statistics see, e.g., Refs. [14–17].

Suppose the outcome of an experiment is characterized by a vector of data $\mathbf{x}$, whose probability distribution depends on an unknown parameter (or parameters) $\theta$ that we wish to determine. In Bayesian statistics, all knowledge about $\theta$ is summarized by

the posterior p.d.f. $p(\theta|\mathbf{x})$, whose integral over any given region gives the degree of belief for $\theta$ to take on values in that region, given the data $\mathbf{x}$. It is obtained by using Bayes' theorem,

$$p(\theta|\mathbf{x}) = \frac{P(\mathbf{x}|\theta)\pi(\theta)}{\int P(\mathbf{x}|\theta')\pi(\theta') d\theta'}, \quad (40.37)$$

where $P(\mathbf{x}|\theta)$ is the likelihood function, i.e., the joint p.d.f. for the data viewed as a function of $\theta$, evaluated with the data actually obtained in the experiment, and $\pi(\theta)$ is the prior p.d.f. for $\theta$. Note that the denominator in Eq. (40.37) serves to normalize the posterior p.d.f. to unity.

As it can be difficult to report the full posterior p.d.f. $p(\theta|\mathbf{x})$, one would usually summarize it with statistics such as the mean, mode or median value and covariance matrix. In addition one may construct intervals with a given probability content, as is discussed in Sec. 40.4.1 on Bayesian interval estimation.

#### 40.2.6.1 Priors

Bayesian statistics supplies no unique rule for determining the prior $\pi(\theta)$; this reflects the analyst's subjective degree of belief (or state of knowledge) about $\theta$ before the measurement was carried out. For the result to be of value to the broader community, whose members may not share these beliefs, it is important to carry out a sensitivity analysis, that is, to show how the result changes under a reasonable variation of the prior probabilities.

One might like to construct $\pi(\theta)$ to represent complete ignorance about the parameters by setting it equal to a constant. A problem here is that if the prior p.d.f. is flat in $\theta$, then it is not flat for a nonlinear function of $\theta$, and so a different parametrization of the problem would lead in general to a non-equivalent posterior p.d.f.

For the special case of a constant prior, one can see from Bayes' theorem (40.37) that the posterior is proportional to the likelihood, and therefore the mode (peak position) of the posterior is equal to the maximum-likelihood estimator. The posterior mode, however, will change in general upon a transformation of parameter. One may use as the Bayesian estimator a summary statistic other than the mode, such as the median, which is invariant under parameter transformation. But this will not in general coincide with the MLE.

The difficult and subjective nature of encoding personal knowledge into priors has led to what is called objective Bayesian statistics, where prior probabilities are based not on an actual degree of belief but rather derived from formal rules. These give, for example, priors which are invariant under a transformation of parameters, or ones which result in a maximum gain in information for a given set of measurements. For an extensive review see, e.g., Ref. [18].

Objective priors do not in general reflect degree of belief, but they could in some cases be taken as possible, although perhaps extreme, subjective priors. The posterior probabilities as well therefore do not necessarily reflect a degree of belief. However one may regard investigating a variety of objective priors to be an important part of the sensitivity analysis. Furthermore, use of objective priors with Bayes' theorem can be viewed as a recipe for producing estimators or intervals which have desirable frequentist properties.

An important procedure for deriving objective priors is due to Jeffreys. According to Jeffreys' rule one takes the prior as

$$\pi(\theta) \propto \sqrt{\det(I(\theta))}, \quad (40.38)$$

where

$$I_{ij}(\theta) = -E \left[ \frac{\partial^2 \ln P(\mathbf{x}|\theta)}{\partial \theta_i \partial \theta_j} \right] \quad (40.39)$$

is the Fisher information matrix. One can show that the Jeffreys prior leads to inference that is invariant under a transformation of parameters. One should note that the Jeffreys prior does not in general correspond to one's degree of belief about the value of a parameter. As examples, the Jeffreys prior for the mean $\mu$ of a Gaussian distribution is a constant, and for the mean of a Poisson distribution one finds $\pi(\mu) \propto 1/\sqrt{\mu}$.
Neither the constant nor $1/\sqrt{\mu}$ priors can be normalized to unit area and are therefore said to be *improper*. This can be allowed because the prior always appears multiplied by the likelihood function, and if the likelihood falls to zero sufficiently quickly then one may have a normalizable posterior density.

An important type of objective prior is the reference prior due to Bernardo and Berger [19]. To find the reference prior for a given problem one considers the Kullback-Leibler divergence $D_n[\pi, p]$ of the posterior $p(\theta|x)$ relative to a prior $\pi(\theta)$, obtained from a set of i.i.d. data $x = (x_1, \dots, x_n)$:

$$D_n[\pi, p] = \int p(\theta|x) \ln \frac{p(\theta|x)}{\pi(\theta)} d\theta. \quad (40.40)$$

This is effectively a measure of the gain in information provided by the data. The reference prior is chosen so that the expectation value of this information gain is maximized for the limiting case of $n \to \infty$, where the expectation is computed with respect to the marginal distribution of the data,

$$p(x) = \int p(x|\theta)\pi(\theta) d\theta. \quad (40.41)$$

For a single, continuous parameter the reference prior is usually identical to the Jeffreys prior. In the multiparameter case an iterative algorithm exists, which requires sorting the parameters by order of inferential importance. Often the result does not depend on this order, but when it does, this can be part of a sensitivity analysis. Further discussion and applications to particle physics problems can be found in Ref. [20].

### 40.2.6.2 Bayesian treatment of nuisance parameters

As discussed in Sec. 40.2.2, a model may depend on parameters of interest $\theta$ as well as on nuisance parameters $\nu$, which must be included for an accurate description of the data. Knowledge about the values of $\nu$ may be supplied by control measurements, theoretical insights, physical constraints, etc. Suppose, for example, one has data $y$ from a control measurement which is characterized by a probability $p(y|\nu)$. Suppose further that before carrying out the control measurement one's state of knowledge about $\nu$ is described by an initial prior $\pi_0(\nu)$, which in practice is often taken to be a constant or in any case very broad. By using Bayes' theorem (40.1) one obtains the updated prior $\pi(\nu)$ (i.e., now $\pi(\nu) = \pi(\nu|y)$, the probability for $\nu$ given $y$),

$$\pi(\nu|y) \propto p(y|\nu)\pi_0(\nu). \quad (40.42)$$

In the absence of a model for $P(y|\nu)$ one may make some reasonable but *ad hoc* choices in order to approximate $\pi(\nu)$. For a single nuisance parameter $\nu$, for example, one might characterize the uncertainty by a p.d.f. $\pi(\nu)$ centered about its nominal value with a certain standard deviation $\sigma_\nu$. Often a Gaussian p.d.f. provides a reasonable model for one's degree of belief about a nuisance parameter; in other cases, more complicated shapes may be appropriate. If, for example, the parameter represents a non-negative quantity then a log-normal or gamma p.d.f. can be a more natural choice than a Gaussian truncated at zero. Note also that truncation of the prior of a nuisance parameter $\nu$ at zero will in general make $\pi(\nu)$ nonzero at $\nu = 0$, which can lead to an unnormalizable posterior for a parameter of interest that appears multiplied by $\nu$.

The likelihood function, prior, and posterior p.d.f.s all depend on both $\theta$ and $\nu$, and are related by Bayes' theorem, as usual. Note that the likelihood here only refers to the primary measurement $x$. Once any control measurements $y$ are used to find the updated prior $\pi(\nu)$ for the nuisance parameters, this information is fully encapsulated in $\pi(\nu)$ and the control measurements do not appear further.

One can obtain the posterior p.d.f. for $\theta$ alone by integrating over the nuisance parameters, i.e.,

$$p(\theta|x) = \int p(\theta, \nu|x) d\nu. \quad (40.43)$$

Such integrals can often not be carried out in closed form, and if the number of nuisance parameters is large, then they can be difficult to compute with standard Monte Carlo methods. Markov Chain Monte Carlo (MCMC) techniques are often used for computing integrals of this type (see Sec. 42.6).

## 40.3 Statistical tests

In addition to estimating parameters, one often wants to assess the validity of certain statements concerning the data's underlying distribution. Frequentist hypothesis tests, described in Sec. 40.3.1, provide a rule for accepting or rejecting hypotheses depending on the outcome of a measurement. In significance tests, covered in Sec. 40.3.2, one gives the probability to obtain a level of incompatibility with a certain hypothesis that is greater than or equal to the level observed with the actual data. Goodness-of-fit tests that quantify the general level of compatibility between data and a hypothesis are described in Sec. 40.3.3. In the Bayesian approach, the corresponding procedure is based fundamentally on the posterior probabilities of the competing hypotheses. In Sec. 40.3.4 we describe a related construct called the Bayes factor, which can be used to quantify the degree to which the data prefer one or another hypothesis.

### 40.3.1 Hypothesis tests

A frequentist test of a hypothesis (often called the null hypothesis, $H_0$) is a rule that states for which data values $x$ the hypothesis is rejected. A region of $x$-space called the critical region, $w$, is specified such that there is no more than a given probability under $H_0$, $\alpha$, called the size or significance level of the test, to find $x \in w$. If the data are discrete, it may not be possible to find a critical region with exact probability content $\alpha$, and thus we require $P(x \in w|H_0) \le \alpha$. If the data are observed in the critical region, $H_0$ is rejected.

The data $x$ used to construct a test could be, for example, a set of values that characterizes an individual event. In this case the test corresponds to classification as, e.g., signal or background. Alternatively the data could represent a set of values from a collection of events. Often one is interested in knowing whether all of the events are of a certain type (background), or whether the sample contains at least some events of a new type (signal). Here the background-only hypothesis plays the role of $H_0$, and in the alternative $H_1$ both signal and background are present. Rejecting $H_0$ is, from the standpoint of frequentist statistics, the required step to establish discovery of the signal process.

The critical region is not unique. Its choice should take into account the probabilities for the data predicted by some alternative hypothesis (or set of alternatives) $H_1$. Rejecting $H_0$ if it is true is called a *type-I error*, and occurs by construction with probability no greater than $\alpha$. Not rejecting $H_0$ if an alternative $H_1$ is true is called a *type-II error*, and for a given test this will have a certain probability $\beta = P(x \notin w|H_1)$. The quantity $1 - \beta$ is called the power of the test of $H_0$ with respect to the alternative $H_1$. A strategy for defining the critical region can therefore be to maximize the power with respect to some alternative (or alternatives) given a fixed size $\alpha$.

To maximize the power of a test of $H_0$ with respect to the alternative $H_1$, the Neyman-Pearson lemma states that the critical region $w$ should be chosen such that for all data values $x$ inside $w$, the likelihood ratio

$$\lambda(x) = \frac{f(x|H_1)}{f(x|H_0)} \quad (40.44)$$

is greater than or equal to a given constant $c_\alpha$, and everywhere outside the critical region one has $\lambda(x) < c_\alpha$, where the value of $c_\alpha$ is determined by the size of the test $\alpha$. The hypotheses $H_0$ and $H_1$ must be simple, i.e., they should not contain undetermined parameters.

It is convenient to define the test using a scalar function of the data $x$ called a *test statistic*, $t(x)$, such that the boundary of the critical region is given by a surface of constant $t(x)$. The Neyman-Pearson lemma is equivalent to the statement that the likelihood ratio (40.44) represents the optimal test statistic. It can be difficult in practice, however, to determine $\lambda(x)$, since this requires knowledge of the joint p.d.f.s $f(x|H_0)$ and $f(x|H_1)$.
Often one does not have explicit formulae for these, but rather Monte Carlo models that allow one to generate instances of $x$ that follow the p.d.f.s.

In the case where the likelihood ratio (40.44) cannot be used explicitly, there exist a variety of other multivariate methods for constructing a test statistic that may approach its performance. These are based on machine-learning algorithms that use samples of training data corresponding to the hypotheses in question, often generated from Monte Carlo models. Further information on Machine Learning can be found in Sec. 41 of this Review.

The multivariate algorithms designed to classify events into signal and background types also form the basis of tests of the hypothesis that a sample of events consists of background only. Such a test can be constructed using the distributions of the test statistic $t(x)$ for event classification obtained from a multivariate algorithm such as a Neural Network output. The distributions $p(t|s)$ and $p(t|b)$ for signal and background events, respectively, are used to construct the likelihood ratio of the signal-plus-background hypothesis relative to that of background only. To the extent that the test statistic $t(x)$ approximates the likelihood ratio (or a monotonic function thereof) for individual events given by (40.44), the resulting test of the background-only hypothesis for the event sample will have maximum power with respect to the signal-plus-background alternative (see Ref. [21]).

### 40.3.2 Tests of significance, p-values

A frequentist significance test allows one to assign a numerical value, called the p-value, that reflects the level of agreement between a hypothesis $H_0$ and the observed data $x$. To carry out a significance test one must specify a subset of the data space, whose boundary includes the observed data point $x_{obs}$, that is deemed to have equal or lesser compatibility with $H_0$ relative to $x_{obs}$. The p-value of $H_0$ is the probability, assuming data distributed according to $H_0$, to find $x$ in the region of equal or lesser compatibility.

For continuous data, one can show that the p-value of $H_0$, here called $p_0$, follows a uniform distribution for data generated under $H_0$. Therefore one has $P(p_0 \le \alpha | H_0) = \alpha$ for any constant $\alpha$ in $[0, 1]$. In this way the region of data space where $p_0 \le \alpha$ gives a critical region $w$ of a test of size $\alpha$, as it satisfies $P(x \in w | H_0) \le \alpha$. Thus rejecting $H_0$ if $p_0 \le \alpha$ is equivalent to the hypothesis test of Sec. 40.3.1 above. For the critical region $w$ defined by $p_0 \le \alpha$ one can find the power $P(x \in w | H_1)$ with respect to an alternative hypothesis $H_1$, as defined in Sec. 40.3.1. In general the region of low compatibility with $H_0$ is defined so as to have a high probability for data generated under $H_1$, so that the distribution of $p_0$, assuming data $x$ that follows $H_1$, will be concentrated at low values, and the critical region $w$ derived from $p_0 \le \alpha$ will have a high power with respect to $H_1$.

Similar to construction of a test's critical region in Sec. 40.3.1, finding the p-value is generally done by defining a statistic $t(x)$ whose value corresponds to better or worse levels of agreement with the hypothesis. The hypothesis being tested, $H_0$, will determine the p.d.f. $f(t|H_0)$ for the statistic. If, for example, $t$ is defined such that large values correspond to poor agreement with the hypothesis, then the p-value is

$$p = \int_{t_{\text{obs}}}^{\infty} f(t|H_0) dt, \qquad (40.45)$$

where $t_{\text{obs}} = t(x_{\text{obs}})$ is the value of the statistic obtained in the actual experiment.

Note that the p-value is not the probability of the hypothesis, which in frequentist statistics is not defined. The p-value should also not be confused with the size (significance level) of a test, or the confidence level of a confidence interval (Section 40.4), both of which are pre-specified constants.

When searching for a new phenomenon, one tries to reject the hypothesis $H_0$ that the data are consistent with known (e.g., Standard Model) processes. If the p-value of $H_0$ is sufficiently low, then one is willing to accept that some alternative hypothesis is true. Often one converts the p-value into an effective significance $Z$, defined as an equivalent number of standard deviations of a

Gaussian distributed random variable. In a search for an intrinsically positive signal, i.e., where only upward fluctuations of the estimated rate appear signal like, this is defined as

$$Z = \Phi^{-1}(1 - p). \qquad (40.46)$$

Here $\Phi$ is the cumulative distribution of the standard Gaussian, and $\Phi^{-1}$ is its inverse (quantile) function. In this way, a p-value of $1/2$ gives $Z = 0$, i.e., having an estimated signal rate of zero corresponds to zero significance. If either positive or negative data fluctuations would indicate evidence of the signal, then one defines

$$Z = \Phi^{-1}(1 - p/2). \qquad (40.47)$$

In this case an estimated signal rate of zero gives $p = 1$ and $Z = 0$.

Often in particle physics the level of significance where an effect is said to qualify as a discovery is $Z = 5$, i.e., a $5\sigma$ effect, corresponding to a p-value of $2.87 \times 10^{-7}$. One's actual degree of belief that a new process is present, however, will depend in general on other factors as well, such as the plausibility of the new signal hypothesis and the degree to which it can describe the data, one's confidence in the model that led to the observed p-value, and possible corrections for multiple observations out of which one focuses on the smallest p-value obtained (the “look-elsewhere effect”, discussed in Section 40.3.2.2).

#### 40.3.2.1 Frequentist treatment of nuisance parameters and asymptotic methods for tests

Suppose one wants to test hypothetical values of parameters $\theta$, but the model also contains nuisance parameters $\nu$. To find a p-value for $\theta$ we can construct a test statistic $q_\theta$ such that larger values constitute increasing incompatibility between the data and the hypothesis. Then for an observed value of the statistic $q_{\theta, \text{obs}}$, the p-value of $\theta$ is

$$p_{\theta}(\nu) = \int_{q_{\theta, \text{obs}}}^{\infty} f(q_{\theta}|\theta, \nu) dq_{\theta}, \qquad (40.48)$$

which depends in general on the nuisance parameters $\nu$. In the strict frequentist approach, $\theta$ is rejected only if the p-value is less than $\alpha$ for all possible values of the nuisance parameters.

The difficulty described above is effectively solved if we can define the test statistic $q_\theta$ in such a way that its distribution $f(q_\theta|\theta)$ is independent of the nuisance parameters. Although exact independence is only found in special cases, it can be achieved approximately by use of the profile likelihood ratio. This is given by the profile likelihood from Eq. (40.18) divided by the value of the likelihood at its maximum, i.e., when evaluated with the maximum-likelihood estimators $\hat{\theta}$ and $\hat{\nu}$:

$$\lambda_p(\theta) = \frac{L(\theta, \hat{\nu}(\theta))}{L(\hat{\theta}, \hat{\nu})}. \qquad (40.49)$$

Wilks' theorem [10] states that, providing certain general conditions are satisfied, the distribution of $-2 \ln \lambda_p(\theta)$, under assumption of $\theta$, approaches a $\chi^2$ distribution in the limit where the data sample is very large, independent of the values of the nuisance parameters $\nu$. Here the number of degrees of freedom is equal to the number of components of $\theta$. More details on use of the profile likelihood are given in Refs. [22, 23] and in contributions to the PHYSTAT conferences [24]; explicit formulae for special cases can be found in Ref. [25]. Further discussion on how to incorporate systematic uncertainties into p-values can be found in Ref. [26].

Even with use of the profile likelihood ratio, for a finite data sample the p-value of hypothesized parameters $\theta$ will retain in general some dependence on the nuisance parameters $\nu$. Ideally one would find the the maximum of $p_\theta(\nu)$ from Eq. (40.48) explicitly, but that is often impractical. An approximate and computationally feasible technique is to use $p_\theta(\hat{\nu}(\theta))$, where $\hat{\nu}(\theta)$ are the profiled values of the nuisance parameters as defined in Section 40.2.2.2. The resulting p-value is correct if the true values of the nuisance parameters are equal to the profiled values used; otherwise it could be either too high or too low. This is discussed further in Section 40.4.2 on confidence intervals.