# Laboratory limits on the annihilation or decay of dark matter particles

Teresa Marrodán Undagoitia<EMAIL_ADDRESS>Max-Planck-Institut
für Kernphysik, Postfach 103980, D-69029 Heidelberg, Germany  Werner
Rodejohann<EMAIL_ADDRESS>Max-Planck-Institut für
Kernphysik, Postfach 103980, D-69029 Heidelberg, Germany  Tim Wolf
<EMAIL_ADDRESS>Max-Planck-Institut für Kernphysik, Postfach 103980,
D-69029 Heidelberg, Germany  Carlos E. Yaguna<EMAIL_ADDRESS>Escuela de Física, Universidad Pedagógica y Tecnológica de Colombia, Avenida
Central del Norte 39-115, Tunja, Colombia

###### Abstract

Constraints on the indirect detection of dark matter are usually obtained from
observations of astrophysical objects – the Galactic Center, dwarf galaxies,
M31, etc. Here we propose instead to look for the annihilation or decay of
dark matter particles taking place inside detectors searching _directly_ for
dark matter or in large neutrino experiments. We show that the data from
XENON1T and Borexino set limits on the annihilation and decay rate of dark
matter particles with masses in the keV to few MeV range. All relevant final
states are considered: annihilation into $\gamma\gamma$ and $e^{-}e^{+}$, and
decays into $\gamma\gamma$, $\gamma\nu$ and $e^{-}e^{+}$. The expected
sensitivities in XENONnT, DARWIN, JUNO and THEIA are also computed. Though
weaker than current astrophysical bounds, the laboratory limits (and
projections) obtained are free from the usual astrophysical uncertainties
associated with $J$-factors and unknown backgrounds, and may thus offer a
complementary probe of the dark matter properties. We point out that current
and future (astro)particle physics detectors might also be used to set
analogous limits for different decays and dark matter masses above a few MeV.

##  I Introduction

Constraints on dark matter (DM) properties can be obtained in different ways.
The main methods are direct detection Marrodán Undagoitia and Rauch (2016);
Billard  _et al._ (2021), in which the cross section of DM particles with
nuclei or electrons is measured, and indirect detection Klasen  _et al._
(2015), in which the decay or annihilation rate of DM particles into Standard
Model (SM) particles is measured. An additional method is the production of DM
particles at colliders. In that case the signatures are events with missing
energy, which is attributed to DM particles being produced and escaping the
detector Kahlhoefer (2017).

Direct detection is subject to many sophisticated experiments, using different
target nuclei and measurement techniques. Looking for tiny recoils stemming
from a DM-nucleus or DM-electron scattering, the constrained parameters are
the mass of the DM particle and its cross section with the target particle.
Indirect detection is typically performed in a parasitic way in telescopes
aiming to perform astrophysical measurements. In this case, the DM
annihilation or decay rate in any SM particle is constrained as a function of
the DM mass, from looking directly for those SM particles, photons being
radiated off them, or their decay products. Inherent to such measurements are
astrophysical uncertainties such as the DM density along the line of sight
(called $J$-factor), or SM processes in the observed complex objects that
generate very similar signatures. Objects that are particularly popular in
what regards observing them are the galactic center, dwarf galaxies, or the
Andromeda galaxy.

In this paper, we suggest to search for events where the DM decay or
annihilation takes place inside Earth-based detectors. The requirements for
such a detector are: the ability to reconstruct the energy deposited by
particles (spectroscopy), a large target volume and ideally a very low
background. Our study is focused on two-phase liquid xenon experiments and
organic scintillation detectors. These technologies provide extremely clean
environments, as demonstrated by world-leading direct detection results and
background rates by XENON1T Aprile  _et al._ (2018) and by the measurements of
all solar neutrino components in Borexino Agostini  _et al._ (2019). The
energy depositions expected in these detectors from DM interactions or solar
neutrinos are in the keV or MeV regime, and the experiments are optimized to
look for such signals. Therefore, our constraints on DM annihilation and decay
will be in this mass range. While it is expected that such laboratory limits
on indirect detection rates will be much weaker than the ones from
astrophysical observations of whole galaxies, the approach proposed here has
several advantages: $i)$ the experimental conditions are under control so the
backgrounds are known; $ii)$ the signal does not suffer from astrophysical
uncertainties such as the $J$-factor; $iii)$ significant improvements in the
limits can be foreseen from the upcoming generation of experiments.

We consider annihilation processes into $e^{-}e^{+}$ and $\gamma\gamma$, and
decays into $e^{-}e^{+}$, $\gamma\gamma$ and $\gamma\nu$. We shall analyze
public XENON1T data Aprile  _et al._ (2020a, b) and extrapolate the limits to
upcoming XENONnT Aprile  _et al._ (2020c) and DARWIN Aalbers  _et al._ (2016)
data, making different assumptions on the background rates in those
experiments. The sensitivity of the upcoming LZ experiment Akerib  _et al._
(2015) should be very similar to that of XENONnT. We employ as well the public
data from Borexino Agostini  _et al._ (2018) to derive results, and make
projections for JUNO Djurcic  _et al._ (2015) and THEIA Askins  _et al._
(2020) using their predicted background spectra.

This paper is built up as follows: first, in Sec. II we discuss the
annihilation and decay processes. In Sec. III we describe general features of
signal production and the background spectra used for the analysis. In Sec. IV
we derive limits on existing and expected future data for the signatures under
study, namely DM annihilation and DM decay. We describe these physics
processes in detail and explain the procedure of limit setting. A discussion
of our results is presented in Sec. V.

##  II Annihilation and decay processes

We consider DM annihilation and decay processes that may take place inside
terrestrial detectors. Their respective rates (number of events per unit
volume and unit time) are given by

$\displaystyle\mathcal{R}$ $\displaystyle=\frac{\langle\sigma
v\rangle}{2}\,n_{\rm DM}^{2}=\frac{\langle\sigma
v\rangle}{2m_{\chi}^{2}}\,\rho_{0,\rm DM}^{2}\quad\mathrm{(Annihilation)},$
(1) $\displaystyle\mathcal{R}$ $\displaystyle=\Gamma_{\rm DM}\,n_{\rm
DM}=\frac{\Gamma_{\rm DM}}{m_{\chi}}\,\rho_{0,\rm DM}\quad\mathrm{(Decay)},$
(2)

where $\langle\sigma v\rangle$ and $\Gamma_{\rm DM}$ are respectively the
annihilation and decay rates into a given final state, $m_{\chi}$ is the mass
of the DM particle, and $\rho_{0,\rm DM}$ is the local DM density, which we
take to be $0.3$ GeV/cm3 Green (2017). Our goal is then to set limits on
$\langle\sigma v\rangle$ and $\Gamma_{\rm DM}$ as a function of $m_{\chi}$ for
all possible final states. Obviously, the smaller the DM mass, the larger is
the rate and the stronger the constraint. Given the relatively small volume of
typical detectors, meaningful bounds can be derived only for light DM
particles (keV to MeV). Notice that, besides the particle physics quantities
we are interested in, the event rates depend only on the _local_ DM density.
Therefore, they are insensitive to the usual large uncertainties associated
with the distribution of DM inside astrophysical objects, typically encoded in
the so-called $J$-factors Strigari (2013). Moreover, unlike direct detection
rates, equations (1) and (2) do not depend on the DM velocity distribution.

The experiments considered in this work are sensitive to energy depositions
between a few keV and a few MeV. For DM annihilations and decays, the total
energy deposition is of the same order as $m_{\chi}$, which implies that the
range of DM masses where bounds can be derived is $\mathrm{keV}\lesssim
m_{\chi}\lesssim\mathrm{MeV}$.

The possible final states from annihilation or decay are quite limited for DM
masses below a few MeV, as the only accessible SM particles are photons,
neutrinos, and electrons/positrons. In this analysis, we consider two DM
annihilation reactions (into $\gamma\gamma$ and $e^{-}e^{+}$) and three decay
modes ($\gamma\gamma$, $\gamma\nu$, and $e^{-}e^{+}$). We remain agnostic
about the nature of the DM particle, but note that decays into a $\gamma\nu$
($\gamma\gamma$) final state resemble a keV-scale sterile neutrino Abazajian
(2017); Boyarsky  _et al._ (2019), axion Chadha-Day  _et al._ (2021) or axion-
like particle Irastorza and Redondo (2018) signals. Our results apply to many
of the countless DM candidates proposed so far.

##  III Detector response and background spectra

Charged particles deposit their energy in liquid xenon producing thereby
scintillation and ionization. In dual-phase time-projection chambers, free
electrons can be collected and amplified via proportional scintillation in the
gas phase on top of the liquid target Aprile  _et al._ (2017). Therefore, both
scintillation photons and the ionization signals are read out using photo-
sensors. These two signals can be combined Aprile  _et al._ (2020b) to improve
the energy resolution of the detector. In organic liquid scintillators,
charged particles deposit their energy similarly exciting the medium. In the
de-excitation process, emitted photons propagate towards the walls of the
detector where photo-sensors are located. While in liquid xenon detectors, the
charge signal can be used to identify multiple scatterings, in organic
scintillators this is not possible. Interactions by mono-energetic
$\gamma$-particles appear in the latter as mono-energetic peaks resulting from
their full absorption. In this work, we consider annihilation and decay
reactions in which observable electrons, positrons and $\gamma$-rays are
created. The $\gamma$-rays have to interact in the sensitive volume in order
to leave a measurable signal. For the energy range considered in this work,
they interact either via photo-absorption or via Compton scattering.

Energy resolution:
The finite energy resolution of the detectors is taken into account when
modeling the signals. We parametrize the energy resolution of the
corresponding signals via the equation

$\frac{\sigma(E\,[\textrm{keV}])}{E}=\frac{a}{\sqrt{E\,[\textrm{keV}]}}+b.$
(3)

In XENON1T, the parameters are determined to be $a=31.3$ and $b=0.17$ Aprile
_et al._ (2020b). As our annihilation and decay signatures involve
$\gamma$-rays, we use for those a slightly different parametrization which
takes into account the fact that above $\sim 350$ keV $\gamma$-rays undergo
multiple scatters. To obtain this slightly different parametrization (with the
parameters $a=31$ and $b=0.19$), we fit the single scatter points below 350
keV and the multiple scatters above that energy (data from figure 6 in Ref.
Aprile  _et al._ (2020b)).

For the organic liquid-scintillation detectors considered here (Borexino, JUNO
and THEIA), we parametrize the energy resolution as in equation (3) setting
$b=0$. For the Borexino experiment, we fit the mono-energetic peaks produced
by $\gamma$ calibration sources as shown in figure 17 of Back  _et al._ (2012)
with Gaussian functions. We extract the corresponding resolution $\sigma(E)$
and obtain a fit value for equation (3) of $a=169$. For the upcoming/future
detectors, we take the expected resolution values as found in the literature
corresponding to $3\%/\sqrt{E\,[\textrm{MeV}]}$ and
$7\%/\sqrt{E\,[\textrm{MeV}]}$ for JUNO Djurcic  _et al._ (2015) and THEIA
Askins  _et al._ (2020), respectively. Table 1 summarizes the energy
resolution for all experiments considered here.

Efficiency and exposure:
While for organic scintillators we assume a detector efficiency of 100%, for
the xenon-based detectors we take an energy-dependent efficiency to be
consistent with the published data (figure 2 of Aprile  _et al._ (2020a)). The
latter has an efficiency loss towards low energies (below $\sim 4$ keV) and at
high energies is flat at the level of 90%.

Our analysis derives results based on the 0.65 tonne$\cdot$y data of the
XENON1T experiment Aprile  _et al._ (2020a), and gives expected sensitivities
for its successor XENONnT Aprile  _et al._ (2020c) and the future DARWIN
experiment Aalbers  _et al._ (2016). The signals in these detectors scale with
the exposure, the product of the fiducial mass (mass used for data analysis)
and the measuring time. The considered exposures are 20 tonnes$\cdot$y and 200
tonnes$\cdot$y for XENONnT and DARWIN, respectively. Note however that for DM
annihilation or decay inside these detectors, the volume (and not the mass) is
relevant, and therefore we divide the exposure by the xenon density
($\rho=2.85$ g/cm3). To derive results for Borexino, we take an exposure of
1291 days and 71.3 tonnes corresponding to a total of 252 tonne$\cdot$year
Agostini  _et al._ (2018). Note that for deposited energies above $\sim 3$ MeV
energies, the fiducial volume of these detectors can be enlarged allowing for
significantly larger exposures Agostini  _et al._ (2020a). This is however not
considered in our study. For JUNO and THEIA, we assume 100
ktonne$\,\cdot\,$year and 300 ktonne$\,\cdot\,$year, respectively, in
accordance with the expected exposures mentioned in the corresponding
publications Djurcic  _et al._ (2015); Askins  _et al._ (2020). While the
density for THEIA is assumed to be 1 g/cm3, we take for Borexino 0.89 g/cm3
Agostini  _et al._ (2020b) and for JUNO 0.86 g/cm3, as reported in Djurcic
_et al._ (2015). Table 1 summarizes the parameters employed for each
experiment.

Table 1: Summary of parameters for the experiments considered in this study. Note that for this analysis the exposure given here has to be divided by the detector material’s density.  Experiment |  Energy resolution (at 1 MeV) |  Exposure | Density
---|---|---|---
XENON1T | 1.2% | 0.65 tonne $\cdot$ y | 2.85 g $\cdot\,\text{cm}^{-3}$
XENONnT | 1.2% | 20 tonne $\cdot$ y | 2.85 g $\cdot\,\text{cm}^{-3}$
DARWIN | 1.2% | 200 tonne $\cdot$ y | 2.85 g $\cdot\,\text{cm}^{-3}$
Borexino | 5.3% | 252 tonne $\cdot$ y | 0.89 g $\cdot\,\text{cm}^{-3}$
JUNO | 3% | 100 ktonne $\cdot$ y | 0.86 g $\cdot\,\text{cm}^{-3}$
THEIA | 7% | 300 ktonne $\cdot$ y | 1.0 g $\cdot\,\text{cm}^{-3}$

Background spectra:
In all considered cases, the signal is a mono-energetic peak appearing on top
of the experimental background. We use the published background spectrum from
XENON1T: from Aprile  _et al._ (2020a) at low energies and from Aprile  _et
al._ (2020b) above $200$ keV. The latter data (given in arbitrary units in the
publication) is scaled to fit the low energy rate. We take the single
scattering data of Aprile  _et al._ (2020b) below 350 keV and the spectrum of
multiple interactions above this energy since the probability of the
$\gamma$-ray to undergo multiple scatters increases with energy. This is a
conservative choice since the background for multiple scatterings is higher
than the single scatter one above 350 keV. Figure 1 left shows, as an example,
the signals from the decay $\chi\to\nu+\gamma$, for a 100, 200, and 400 keV DM
particle (blue, green, and red, respectively) on top of the background
spectrum of XENON1T (in gray). Examples for the same decay channel are shown
in figure 1 right for DM masses corresponding to 1.5 MeV, 3.3 MeV, and 5.0
MeV. The underlying DM decay rate for all examples is
$\Gamma=10^{-13}\,\textrm{s}^{-1}$.

Figure 1: The XENON1T background (from Aprile  _et al._ (2020b, a) as
described in the text) and some typical signals. The left panel displays the
background in log-scale between $1$ keV and $3$ MeV, while the right panel
shows the same background in linear scale and only between $400$ keV and $3$
MeV. Some representative signals from a DM decay via $\chi\to\nu+\gamma$ are
shown, on top of the background, for different DM masses.

The sensitivity of the experiment depends on the specific DM mass due to the
features of the background spectrum and the expected signal rates (see
equations (1, 2)). Several lines are visible in figure 1 originating from
backgrounds in the detector, for example the 41 keV line from 83mKr or the
1.17 MeV and 1.33 MeV lines from 60Co decays. For masses corresponding to
peaks in the background spectrum, the sensitivity can be up to an order of
magnitude lower than for masses in-between.

To calculate the sensitivity of future DM detectors to DM annihilation and
decay, we need to make some assumptions on their background levels. For
XENONnT, we take the background shape from XENON1T and assume a background
suppression factor of 8. We study however also the effect of a pessimistic
(factor 5) and more optimistic (factor 10) background suppression. Similarly,
we employ a factor 80 as suppression factor for the DARWIN background compared
to the background of XENON1T and test the effect of factors 50 and 100.

We take as well the published background spectrum of Borexino from Agostini
_et al._ (2018) to derive results. The left panel of figure 2 shows this
spectrum together with some representative decay lines ($\chi\to\nu+\gamma$)
for a DM mass of 1.5, 3.0, and 5.0 MeV (in blue, green, and red respectively)
and the decay rate is either $10^{-14}\,$s-1 (blue, green) or $10^{-15}\,$s-1
(red). The structure of the spectrum is given by the sum of the various solar
neutrino components (7Be and 8B for instance) and some backgrounds from
natural radioactivity in the detector materials. The sensitivity for a signal
depositing around $\sim 500$ keV is significantly affected by the 210Po
$\alpha$-decay peak at this energy. The energy threshold of Borexino is
significantly higher than for the xenon-based experiments because the
background from 14C present in the organic liquid becomes dominant below $\sim
150$ keV (see the rise of the rate in figure 2).

Figure 2: Left: Borexino background spectrum (from Agostini  _et al._ (2018))
including some typical signals. The right panel shows the comparison of the
Borexino spectrum to the expected spectra of JUNO and THEIA from Djurcic  _et
al._ (2015) and Askins  _et al._ (2020), respectively.

For JUNO, we use the published expected spectrum (top panel of figure 1.10 in
Djurcic  _et al._ (2015)) which starts close to the end-point of the 14C
spectrum and covers up to 1.8 MeV where the background from 11C dominates. In
the case of THEIA, we take the expected background spectrum for a 25 ktonne
target from figure 8 in Askins  _et al._ (2020). This spectrum starts at a bit
higher energy and ranges up to 6 MeV. The peak at $\sim 4$ MeV originates from
208Tl of the thorium chain inside the scintillator including the 2.6 MeV
$\gamma$-line and the energy deposition from the associated beta decay. Below
this energy, signals from 7Be neutrinos and 40K background decays determine
the shape of the spectrum. The used spectra for JUNO and THEIA are shown in
the right panel of figure 2.

The precise knowledge of the background is one of the main advantages of the
laboratory searches we are proposing here. In astrophysical environments, it
is essentially impossible to measure the background. Thus, given a putative
signal, it is extremely difficult to ascertain that it is actually due to DM.
The long-standing controversies regarding the Galactic Center excess (see e.g.
Goodenough and Hooper (2009); Bartels  _et al._ (2016)) and the 3.5 keV line
(see e.g. Boyarsky  _et al._ (2015); Jeltema and Profumo (2015)), for
instance, are mostly due to an uncertain background.

##  IV Results

In this section, we show the obtained limits for the annihilation cross
section or the decay rate for each of the final states considered, and the
projected sensitivity of future experiments. A $\chi^{2}$-method is employed
to test the signal hypothesis. For this purpose, we take the background
spectrum as described in section III and test for each DM mass the allowed
signal-strength (annihilation cross section or decay rate) of the signal
models discussed in section II. In all cases, the signal is modeled as a
Gaussian peak with a width given by the corresponding energy resolution.
Details regarding the features of the DM annihilation and decay signals can be
found in section IV.1 and section IV.2. We employ a $\chi^{2}$ probability
density function (pdf) with one degree of freedom,

$p\left(x\right)=\frac{1}{\sqrt{2\pi}}~{}x^{-1/2}\cdot e^{-x/2},$ (4)

to approximate the log-likelihood ratio test-statistics of a binned Poisson
hypothesis test with one free parameter in the signal model Cowan  _et al._
(2011).

As a threshold for the upper limit, we take the 90% quantile of the $\chi^{2}$
pdf in equation (4). This threshold is independent of the tested signal
strength (Wilk’s theorem) which was verified for the obtained limits. The
threshold is compared to the computed $\chi^{2}$ of the signal-hypothesis
under consideration,

$\chi^{2}=\sum_{i=0}^{N_{\text{bins}}}\frac{(O_{i}-E_{i})^{2}}{E_{i}},$ (5)

where $O_{i}$ corresponds to the bin-content when a signal is present and
$E_{i}$ corresponds to the bin-content without signal with a 1 keV binning.
This approximates a log-likelihood ratio scan of a binned Poisson hypothesis
test and takes into account the statistical uncertainty in each bin. The
intersection between the 90% quantile of the $\chi^{2}$ pdf and the scan
through $\chi^{2}$-values per mass point as computed in equation (5) denotes
the upper limit that we set.

###  IV.1 Dark Matter Annihilation

$\boldsymbol{\chi+\chi\to\gamma+\gamma}$: The $\gamma$-rays produced in this
annihilation reaction can only be observed indirectly once they interact with
electrons in the medium. We therefore observe the electrons produced via the
photo-effect or Compton effect. Figure 3 shows schematically a possible event
topology.

Figure 3: Schematic of the annihilation reaction $\chi+\chi\to\gamma+\gamma$,
where one of the $\gamma$’s (bottom) undergoes photo-absorption and the second
(top) experiences a Compton scattering before being absorbed. The electrons
created deposit their kinetic energy leading to the detected signals.

The mass of the DM particles corresponds to the energy of the $\gamma$-rays
and therefore the total energy deposited is:

$E_{{\rm dep}}^{\chi\chi\to\gamma\gamma}=2\cdot m_{\chi}.$ (6)

For low energies, an electron is emitted mostly in the photo-effect process
and deposits its energy after traveling a short distance. With increasing
energy, the probability to undergo multiple scatterings grow. Indeed, the
XENON1T data shows (figure 6 in Aprile  _et al._ (2020b)) that the background
rate of multiple scatterings is higher for energies above 350 keV.
Consequently for liquid xenon detectors, the single scattering spectrum is
employed below 350 keV and the multiple scattering spectrum above that energy,
as explained in Sec. III.

The results for this case are shown in the left panel of figure 4. The region
above the solid red line is excluded by the data from XENON1T. The dashed
(orange) and dotted (magenta) lines show instead the expected sensitivities of
XENONnT and DARWIN which are predicted to strengthen XENON1T limit by factors
of about $20$ and $200$ respectively over the entire mass range.

Figure 4: Limits and prospects for $\chi+\chi\rightarrow\gamma+\gamma$ (left)
and $\chi+\chi\rightarrow e^{-}+e^{+}$ (right). The red line (solid) shows the
derived limit from XENON1T whereas the orange (dashed) and magenta (dotted)
lines correspond to the expected sensitivities of XENONnT and DARWIN
respectively. The computed limits for Borexino, THEIA, and JUNO are shown in
aquamarine (solid), green (dashed-dotted), and blue (dashed) respectively.
Notice that in the right panel the range of DM masses is smaller (because
$m_{\chi}>m_{e}$) and the scale is linear rather than logarithmic. For
comparison, the lower panels show the corresponding bounds from M31 Ng  _et
al._ (2019) or the CMB Slatyer (2016).

As expected from equation (1), the limits decrease quadratically with the DM
mass. On top of this general trend, the effect of the bumps present in the
background is clearly visible. Current limits of XENON1T on $\langle\sigma
v\rangle$, then, lie between $10^{-21}\,\mathrm{cm^{3}\,s^{-1}}$ at a mass of
$4$ keV and $10^{-15}\,\mathrm{cm^{3}\,s^{-1}}$ at a $2$ MeV mass.

The limits derived from Borexino data, significantly more stringent, are shown
by a continuous aquamarine line in figure 4. Borexino constrains
$\langle\sigma v\rangle$ to be below between
$10^{-19}\,\mathrm{cm^{3}\,s^{-1}}$ at a DM mass of 1 MeV. The expected limits
of JUNO and THEIA are shown in blue and green dashed-lines, respectively, and
are expected to be one to two orders of magnitude stronger mainly due to their
higher exposures. The shape of the limit curves for these experiments can also
be understood by the shape of the background spectra in figure 2 and the
dependence on $1/m_{\chi}^{2}$ of the annihilation rate. The larger volume of
Borexino, THEIA and JUNO (see table 1) implies that the limits are stronger,
but note that the xenon-based experiments can probe much smaller DM masses,
because of their much smaller energy threshold. The liquid xenon experiments
need a $\sim$ keV energy threshold in order to directly detect DM. Above $\sim
100$ keV and depending on the DM mass, the constraints from Borexino can be
$(3-5)$ orders of magnitude better than XENON1T.

$\boldsymbol{\chi+\chi\to e^{-}+e^{+}}$: For the annihilation into electron
and positron, the two charged leptons deposit their energy in the medium. In
addition, once the $e^{+}$ is slowed down, it annihilates with an electron
emitting two 511 keV $\gamma$-rays. Figure 5 shows a scheme of the
annihilation reaction.

Figure 5: Schematic of the annihilation reaction $\chi+\chi\to e^{-}+e^{+}$.
The positron deposits its energy and annihilates with an electron of the
medium producing two 511 keV $\gamma$-rays. These interact promptly with the
medium depositing also their energy.

In this case, the minimum DM mass that can be tested is 511 keV as this is the
mass required to create an electron-positron pair. The total energy deposited
corresponds also to

$E_{{\rm dep}}^{\chi\chi\to e^{-}e^{+}}=2\cdot m_{\chi},$ (7)

because the energy required for the creation of the two leptons is obtained
back when the positron annihilates with an electron at rest. For the liquid
xenon case, we employed a simplified simulation to study the energy resolution
resulting from the depositions of the $e^{-}e^{+}$ and the $\gamma$’s. We
found that the resolution can be modeled by taking the multiple scattering
curve for the total energy as described in Sec. III. With this choice, we
obtain conservative results as the actual resolution is slightly better.

The right panel of figure 4 displays our results for this case. Notice that
the mass range is much smaller and that the scale is linear rather than
logarithmic. Since there is a kinematic threshold for the DM mass in this
channel, Borexino, JUNO and THEIA give better limits than the xenon-based
experiments for all masses. The fact that the limits turn out to be
essentially flat over the entire mass range is a consequence of the decreasing
background – see the right panel of figures 1 and 2. The limits and
projections obtained are seen to be comparable to those from the
$\gamma\gamma$ final state at the same masses.

###  IV.2 Dark Matter Decay

We have seen in the previous subsection how the shape of the limit curves can
be understood from the shape of the background spectra, the volume of the
experiments and the $1/m_{\chi}^{2}$ mass-dependence of the rate. For decay,
everything can be interpreted in analogy, with the exception that the rate
scales with $1/m_{\chi}$. For DM decay, there is an unavoidable limit in terms
of the lifetime $t_{0}$ of the Universe. Some of our limits will lie above
this value and are hence unphysical. This may be expected giving the
comparably small volume of the terrestrial experiments. Our approach is to use
existing and expected data to set limits, so we nevertheless will display the
results above this unavoidable value.

$\boldsymbol{\chi\to\nu+\gamma}$: This is the most important decay mode for
sterile neutrinos, the archetype DM candidate at keV masses. In this decay
mode, the neutrino escapes the detector and only the $\gamma$-ray is detected.
The energy deposited by the $\gamma$-ray corresponds in this case to half of
the DM mass:

$E_{{\rm dep}}^{\chi\to\nu\gamma}=m_{\chi}/2.$ (8)

Our results for this case are shown in the top-left panel of figure 6. The
limits on $\Gamma_{\rm DM}$ tend to weaken linearly with the DM mass,
according to equation (2) and up to background effects. The actual limit on
$\Gamma_{\rm DM}$ from XENON1T varies from about $10^{-17}\,$s-1 at a mass of
few keV to $10^{-13}\,$s-1 at a mass of few MeV. In DARWIN, these limits may
be improved by about three orders of magnitude. Also here, as the xenon-based
experiments have a lower energy threshold, they can test much lower DM masses
while Borexino, JUNO and THEIA provide in turn better limits for the higher
energy range. For instance, Borexino gives a limit of $10^{-19}\,$s-1 at a
mass of 1 MeV, which can be improved by an order of magnitude by THEIA and
JUNO.

Figure 6: Limits and prospects for $\chi\rightarrow\nu+\gamma$ (top left),
$\chi\to\gamma+\gamma$ (top right) and $\chi\to e^{-}+e^{+}$ (bottom). The
solid lines show the results from XENON1T (red) and Borexino (aquamarine). The
projections for XENONnT and DARWIN are shown in dashed (orange) and dotted
(magenta) lines respectively, while those for THEIA and JUNO are displayed as
dashed-dotted (green) and dashed (blue) lines. Notice that in the bottom panel
the range of DM masses is smaller (because $m_{\chi}>2m_{e}$) and the scale is
linear rather than logarithmic. The horizontal dotted (gray) line corresponds
to a lifetime equal to the age of the Universe, $t_{0}$. For comparison, the
lower panels show the respective bounds from M31 Ng  _et al._ (2019) or the
CMB Slatyer and Wu (2017)

.

$\boldsymbol{\chi\to\gamma+\gamma}$: The mass of the DM particle is
transferred to kinetic energy of the $\gamma$-rays. Therefore, the total
energy deposition of the $\gamma$’s corresponds to the DM mass:

$E_{{\rm dep}}^{\chi\to\gamma\gamma}=m_{\chi}.$ (9)

The top-right panel of figure 6 shows the current limit and future prospects
for this decay. The results compare to each other similarly to the results of
the $\gamma+\nu$ decay channel.

$\boldsymbol{\chi\to e^{-}+e^{+}}$: This decay mode is similar to the
annihilation into $e^{-}e^{+}$. The two leptons deposit their kinetic energy
and once stopped, the positron annihilates with an electron producing two 511
keV $\gamma$-rays. The total energy deposition corresponds in this case to the
DM mass:

$E_{{\rm dep}}^{\chi\to e^{-}e^{+}}=m_{\chi}.$ (10)

While the limit of XENON1T on $\Gamma_{\rm DM}$ lies at the level of
$10^{-13}\,$s-1, for Borexino the limit is significantly lower at the level of
a few times $10^{-17}\,$s-1. The future JUNO and THEIA detectors would improve
Borexino results by an order of magnitude.

For all annihilation and decay channels described above, the effect of the
background on the results of the upcoming xenon experiments has been studied.
We assume background suppression factors of 8 and 80 for XENONnT and DARWIN,
respectively. This is also the background level used for all results of the
figures above. We test the effect of different background level into the
results and find a small changes in the obtained limits. For a pessimistic
scenario (factors 5 and 50), the limits worsen by 12%, while for an optimistic
scenario (factors 10 and 100), the limits improve by 25% for both XENONnT and
DARWIN.

##  V Discussion and Conclusion

Large dark-matter and neutrino detectors have a rich experimental programs
beyond their main scientific goals. We have reported here on a new possibility
in which the decay or annihilation of DM particles may take place inside the
Earth-based detectors. Focusing on XENON1T and Borexino as well as the future
XENONnT, JUNO, DARWIN and THEIA, we have set current and prospective future
limits on annihilation/decay rates for DM masses in the keV to few MeV range.
The relevant final states thus involve electrons, positrons and photons.
Usually, cosmological considerations Slatyer (2016) or observations of
astrophysical objects Ng  _et al._ (2019) are used to set such limits, and
they turn out to be more stringent, by several orders of magnitude, than the
ones we derived. Nevertheless, the limits obtained here under well-controlled
laboratory conditions offer the benefits of a well-known background
composition. Moreover, the signal does not suffer from uncertainties
associated with the DM distribution – the $J$\- and $D$-factors. While the
astrophysical and cosmological limits are not expected to improve
significantly, upcoming direct detection and neutrino experiments will improve
the laboratory limits and close the gap considerably. Indeed at low energies
(below $\sim 100$ keV), the upcoming XENONnT and DARWIN experiments will
improve the limit of XENON1T by at least an order a magnitude each. Similarly
above this energy, the larger detectors JUNO and THEIA will improve the
results of Borexino by an order of magnitude.

Similar studies could be performed in argon dark-matter detectors like
DarkSide and Argo Aalseth  _et al._ (2018). Moreover, Super-Kamiokande Fukuda
_et al._ (2003) and Hyper-Kamiokande Abe  _et al._ (2018) have much larger
volumes and could be employed for similar studies at MeV-GeV energies. DUNE
Abi  _et al._ (2020) or large-scale experiments looking for neutrinoless
double beta decay such as KamLAND-Zen Gando  _et al._ (2016) are other up-
coming facilities with features that allow to set limits along the lines
proposed here. The huge detectors at LHC, which are optimized for energies in
the TeV-range, could also provide limits for some final states and DM masses.
All such studies, however, lie beyond the scope of present paper. Our focus
was here on final states involving photons, electrons and positrons. In
principle, for larger DM masses also other particles could be produced, for
instance $\mu^{-}$ and $\mu^{+}$. While being a spectacular signature, energy
deposition by the muons followed by a delayed electron/positron signal from
the two muon decays), the background and efficiencies need to be carefully
evaluated before analyzing this further. The same holds for other final
states. Of course, with increasing DM mass the amount of particles in the
detector volume decreases.

In summary, we have suggested a new method to test DM candidates using direct-
detection and neutrino experiments as indirect detection laboratory
experiments. We obtained present and future constraints on the decay and
annihilation rates of DM particles.

###### Acknowledgements.

We thank our colleague Michael Wurm for his valuable help. CY would like to
thank Kenneth Garcia for his coding help in the initial stages of this work.

## References

  * Marrodán Undagoitia and Rauch (2016) T. Marrodán Undagoitia and L. Rauch, J. Phys. G 43, 013001 (2016), arXiv:1509.08767 .
  * Billard  _et al._ (2021) J. Billard _et al._ , (2021), arXiv:2104.07634 .
  * Klasen  _et al._ (2015) M. Klasen, M. Pohl, and G. Sigl, Prog. Part. Nucl. Phys. 85, 1 (2015), arXiv:1507.03800 .
  * Kahlhoefer (2017) F. Kahlhoefer, Int. J. Mod. Phys. A 32, 1730006 (2017), arXiv:1702.02430 .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 .
  * Agostini  _et al._ (2019) M. Agostini _et al._ (Borexino), Phys. Rev. D 100, 082004 (2019), arXiv:1707.09279 .
  * Aprile  _et al._ (2020a) E. Aprile _et al._ (XENON), Phys. Rev. D 102, 072004 (2020a), arXiv:2006.09721 .
  * Aprile  _et al._ (2020b) E. Aprile _et al._ (XENON), Eur. Phys. J. C 80, 785 (2020b), arXiv:2003.03825 .
  * Aprile  _et al._ (2020c) E. Aprile _et al._ (XENON), JCAP 11, 031 (2020c), arXiv:2007.08796 .
  * Aalbers  _et al._ (2016) J. Aalbers _et al._ (DARWIN), JCAP 1611, 017 (2016), arXiv:1606.07001 .
  * Akerib  _et al._ (2015) D. S. Akerib _et al._ (LZ), (2015), arXiv:1509.02910 .
  * Agostini  _et al._ (2018) M. Agostini _et al._ (Borexino), Nature 562, 505 (2018).
  * Djurcic  _et al._ (2015) Z. Djurcic _et al._ (JUNO), (2015), arXiv:1508.07166 .
  * Askins  _et al._ (2020) M. Askins _et al._ (Theia), Eur. Phys. J. C 80, 416 (2020), arXiv:1911.03501 .
  * Green (2017) A. M. Green, Phys. Rev. D 96, 043020 (2017), arXiv:1705.10818 .
  * Strigari (2013) L. E. Strigari, Phys. Rept. 531, 1 (2013), arXiv:1211.7090 .
  * Abazajian (2017) K. N. Abazajian, Phys. Rept. 711-712, 1 (2017), arXiv:1705.01837 .
  * Boyarsky  _et al._ (2019) A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, Prog. Part. Nucl. Phys. 104, 1 (2019), arXiv:1807.07938 .
  * Chadha-Day  _et al._ (2021) F. Chadha-Day, J. Ellis, and D. J. E. Marsh, (2021), arXiv:2105.01406 .
  * Irastorza and Redondo (2018) I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018), arXiv:1801.08127 .
  * Aprile  _et al._ (2017) E. Aprile _et al._ (XENON), Eur. Phys. J. C 77, 881 (2017), arXiv:1708.07051 .
  * Back  _et al._ (2012) H. Back _et al._ (Borexino), JINST 7, P10018 (2012), arXiv:1207.4816 .
  * Agostini  _et al._ (2020a) M. Agostini _et al._ (Borexino), Phys. Rev. D 101, 062001 (2020a), arXiv:1709.00756 .
  * Agostini  _et al._ (2020b) M. Agostini _et al._ (Borexino), Phys. Rev. D 101, 012009 (2020b), arXiv:1909.02257 .
  * Goodenough and Hooper (2009) L. Goodenough and D. Hooper, (2009), arXiv:0910.2998 .
  * Bartels  _et al._ (2016) R. Bartels, S. Krishnamurthy, and C. Weniger, Phys. Rev. Lett. 116, 051102 (2016), arXiv:1506.05104 .
  * Boyarsky  _et al._ (2015) A. Boyarsky, J. Franse, D. Iakubovskyi, and O. Ruchayskiy, Phys. Rev. Lett. 115, 161301 (2015), arXiv:1408.2503 .
  * Jeltema and Profumo (2015) T. E. Jeltema and S. Profumo, Mon. Not. Roy. Astron. Soc. 450, 2143 (2015), arXiv:1408.1699 .
  * Cowan  _et al._ (2011) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), [Erratum: Eur.Phys.J.C 73, 2501 (2013)], arXiv:1007.1727 .
  * Ng  _et al._ (2019) K. C. Ng, B. M. Roach, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos, and D. R. Wik, Phys. Rev. D 99, 083005 (2019), arXiv:1901.01262 .
  * Slatyer (2016) T. R. Slatyer, Phys. Rev. D 93, 023527 (2016), arXiv:1506.03811 .
  * Slatyer and Wu (2017) T. R. Slatyer and C.-L. Wu, Phys. Rev. D 95, 023010 (2017), arXiv:1610.06933 [astro-ph.CO] .
  * Aalseth  _et al._ (2018) C. E. Aalseth _et al._ (DarkSide-20k), Eur. Phys. J. Plus 133, 131 (2018), arXiv:1707.08145 .
  * Fukuda  _et al._ (2003) Y. Fukuda _et al._ (Super-Kamiokande), Nucl. Instrum. Meth. A 501, 418 (2003).
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 .
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE), JINST 15, T08010 (2020), arXiv:2002.03010 .
  * Gando  _et al._ (2016) A. Gando _et al._ (KamLAND-Zen), Phys. Rev. Lett. 117, 082503 (2016), [Addendum: Phys.Rev.Lett. 117, 109903 (2016)], arXiv:1605.02889 .


Among all of the first-year graduate students enrolled in doctoral-granting
physics departments, the percentage of female and racial minority students has
remained unchanged for the past 20 years. The current graduate program
admissions process can create challenges for achieving diversity goals in
physics. In this paper, we will investigate how the various aspects of a
prospective student’s application to a physics doctoral program affect the
likelihood the applicant will be admitted. Admissions data was collected from
a large, Midwestern public research university that has a decentralized
admissions process and included applicants’ undergraduate GPAs and
institutions, research interests, and GRE scores. Because the collected data
varied in scale, we used supervised machine learning algorithms to create
models that predict who was admitted into the PhD program. We find that using
only the applicant’s undergraduate GPA and physics GRE score, we are able to
predict with 75% accuracy who will be admitted to the program.

# Using machine learning to understand physics graduate school admissions

Nicholas T. Young  Department of Physics and Astronomy, Michigan State
University, East Lansing, Michigan 48824  Department of Computational
Mathematics, Science, and Engineering, Michigan State University, East
Lansing, Michigan 48824  Marcos D. Caballero  Department of Physics and
Astronomy, Michigan State University, East Lansing, Michigan 48824  Center for
Computing in Science Education & Department of Physics, University of Oslo,
N-0316 Oslo, Norway  CREATE for STEM Institute, Michigan State University,
East Lansing, Michigan 48824

##  I Introduction

Despite other science, technology, engineering, and mathematics (STEM) fields
becoming more diverse over the past few decades, physics has lagged behind
with only 20% of bachelor’s degrees awarded to women and only 11% awarded to
racial minorities [1]. These numbers do not improve when considering graduate
degrees where 20% of doctoral degrees are granted to women and 7% are granted
to racial minorities [1]. While this underrepresentation has both enrollment
and retention causes, this paper will focus on the factors that may affect
enrollment in physics graduate programs.

When considering enrollment in physics PhD programs, prior work has found that
minority students in physics are less likely to apply to programs if they feel
that they will not be admitted based on low GPA or GRE scores, or lack of
research experience [2]. Further, given that many graduate programs have
application fees, financial concerns may prevent students from applying to
graduate programs that they believe they will not be admitted to. Therefore,
it is important to understand what matters when applying to physics graduate
programs while acknowledging that many factors that are not easily
quantifiable matter.

Previous research into graduate admissions in physics has tended to take a
broad approach, characterizing the graduate admissions process across the
United States, both for master’s and PhD programs [3, 4] or focusing on a
specific subset of universities such as elite universities [5]. These studies
find that faculty consider numerical measures such as undergraduate GPA and
GRE scores most important in the admissions process and have been conducted by
either observing the admissions process or by surveying faculty about what
they believe to be most important in the admissions process. More recent work
in physics graduate admissions has explored applicant perceptions of the
various components of the application [6]. Missing in this analysis is an
investigation of the actual applications of prospective physics graduate
students. To our knowledge, there has only been one such study [7].

Given that applications to graduate programs consist of numerical data such as
GPA and GRE scores, categorical data such as gender, race, and ethnicity, and
open-ended data such as letters of recommendation and personal statements,
graduate admissions is an ideal target for machine learning. Indeed, machine
learning approaches to understanding graduate admissions have been employed in
computer science to study self-reported admissions data [8] and to streamline
the review process [9]. Machine learning methods have also been employed more
broadly in higher education admissions to predict which admitted students will
accept an offer to attend a small liberal arts school [10, 11] and to predict
which students are likely to be admitted and to complete their MBA [12].

The goal of this work to further the study of graduate admissions in physics
by analyzing the applications using a machine learning approach. In addition,
unlike other studies in physics graduate admissions, this work represents a
case study of a single institution rather than a broad look at the graduate
admissions landscape. However, since physics is regarded as a high consensus
discipline, that is, there is large agreement about what counts as legitimate
admissions practices [13], we expect our results to generalize to similar
doctoral programs. Future work will investigate the robustness of the findings
presented here.

##  II Methods

###  II.1 Data

The data used in this study comes from admissions data recorded from 512
domestic applicants to a physics and astronomy PhD program at a large,
research-intensive, Midwestern public university from 2014-2017. Domestic and
international applicants do not undergo the same review process and hence we
only analyze applications from domestic students here. The admissions process
is unique at this university in that the applications are not reviewed by a
central committee but rather members of the subdiscipline in which the student
expresses interest. For example, if the student intends to study high energy
physics, only high energy physicists would review the application and decide
whether to accept the student. The data include the applicant’s undergraduate
institution and grade point average (GPA), their general and physics GRE
scores, and their physics subdiscipline of interest. As laws of the state
where this university is located prohibit giving preferential treatment to
applicants based on gender, race, ethnicity, or national origin, demographics
were not recorded by the physics and astronomy department and hence, are not
available to us. Overall, 48% of the domestic applicants were offered
admission into the program.

###  II.2 Describing Undergraduate Institutions

Because the name of the undergraduate institution in itself does not provide
useful information to an algorithm, we created new factors to describe
characteristics of the institutions. To describe the overall institution, we
classified each institution as public or private, whether it is a minority
serving institution (MSI), the region of the country it is located in (such as
Northeast, Southwest, etc.), and the Barron’s selectivity of the institution,
which describes how selective the undergraduate program is and serves as a
proxy for prestige. Classifications for the first three categories were taken
from the most recent Carnegie Rankings [14] while the Barron’s classification
came from Barron’s Profiles of American Colleges. Because the overall
reputation of the applicant’s undergraduate university may not describe the
physics program at that university, we also included factors related to the
physics program such as the highest physics degree offered at the university
and the size of the undergraduate program and PhD program if applicable. The
size of the undergraduate and PhD programs were determined by the number of
graduates of the program in a typical year. The programs were then classified
as small, medium-small, medium-large, or large based on which quartile they
fell into. We used the Roster of Physics Departments with Enrollment and
Degree Data to collect this data [15]. All factors appearing in our model are
shown in Table 1 and include the scale of measurement.

Table 1: Variables used in our model and their scale of measurement Factor | Measurement Scale
---|---
Undergraduate GPA | Continuous
Verbal GRE score | Continuous
Quantitative GRE score | Continuous
Written GRE score | Continuous
Physics GRE score | Continuous
Proposed research area | Categorical
Application year | Categorical
Barron’s selectivity | Categorical
Region of applicant’s undergraduate institution | Categorical
Type of physics program at applicant’s undergraduate institution | Categorical
Size of undergraduate physics program at applicant’s undergraduate institution | Categorical
Size of doctoral physics program at applicant’s undergraduate institution | Categorical
Applicant attended a minority serving institution | Binary
Public or Private | Binary
Output variable: admitted status | Binary

###  II.3 Justifying our choices of institutional factors

Prior work has documented university pedigree is often considered in the
application process because institutional quality is assumed to be a proxy for
student quality [5, 16]. Here, we measure institutional quality by Barron’s
selectivity and public or private status, with the assumption that physics
faculty view private universities as more prestigious than public
universities. We include region of the applicant’s undergraduate university to
account for the fact that the institution being studied is a public university
and may therefore show a preference for students from the surrounding region.

Prior work has also found faculty exhibit a tendency to admit students like
themselves, though it is more common among academics who graduated from elite
institutions [5]. Therefore, it is not unreasonable to expect that faculty may
prefer to admit students who followed similar paths as they did, meaning
students from large, doctoral institutions may be more likely to be admitted
than students from smaller institutions. Additionally, we use the size of the
undergraduate and PhD programs as proxies for the perceived prestige of the
physics department, assuming a more prestigious physics department attracts
more students and hence graduates more students.

###  II.4 Random Forest Model

To analyze our data, we used the conditional inference forest algorithm, a
variant of the random forest algorithm shown to be less biased when the data
includes both continuous and categorical variables [17] such as those used in
our model (see table 1). Random forest models in general are ensembles of
individual decision trees, which use binary splits of the input features in
order to make a prediction. The predictions are then averaged over the
individual trees to obtain the overall prediction of the random forest.

While there are multiple metrics used to assess random forest and other
machine learning models, two of the most common are the accuracy and the area
under the curve (AUC). The accuracy is simply the proportion of correct
predictions made by the model. To ensure that the accuracy isn’t inflated by
overtraining, only a fraction of the available data is used to construct the
model while the rest is used to test the predictive power. It is this
remaining data that is used to calculate the accuracy of the model. The AUC is
defined as the area beneath the receiver operator curve of the model, which
visualizes the false positive rate against the true positive rate and varies
between 0.5 and 1, with values greater than 0.7 signifying a good model [18].
This area describes the proportion of positive cases that are ranked above
negative cases in the data set. For example, for our data, the AUC represents
the proportion of all random pairs of admitted and not-admitted applicants in
which the admitted applicant is classified as admitted and the not-admitted
applicant is classified as not-admitted.

In addition to making predictions, the random forest algorithm can determine
the importance of each feature to the model, referred to as the feature
importance. We use the AUC feature importance [19] as it is less biased when
input features differ in scale (as do our factors listed in table 1) and when
the predicted variable is not split evenly between the two outcomes. As the
feature importances are not assumed to follow any statistical distribution,
there is no simple way to apply the idea of statistical significance to
feature importances. We therefore apply the recursive backward elimination
technique described in [20] to determine which features are important and
which are not. We will refer to the features determined to be important via
this process as meaningful features. For more information about random forest
models, biases, and feature importance measures, see Young et al. [21].

We chose to apply a random forest model instead of a more traditional
technique for classifying data such as logistic regression (as used by Attiyeh
and Attiyeh [22] and Posselt et al [7] to study graduate admissions) due to
these feature importances. As feature importances measure all factors on the
same scale, that is how much they change the area under the curve, factors of
otherwise different scales can be compared. This is in contrast to logistic
regression where the odds ratio for a continuous variable would measure the
change in odds for a unit increase in the variable while the odds ratio for a
categorical or binary variable measures the change in odds relative to a
reference group.

We used R [23] for our analysis and the party package [24, 17, 25] to create a
conditional inference forest model. We used 500 trees to build our forest and
70% of our data to train the model. We ran our model 30 times, randomly
selecting 70% of our data for training each time, and averaged the feature
importances over runs so that the resulting distribution of individual feature
importances would be approximately normal. As the conditional inference forest
algorithm has routines built in to handle missing data [26], applicants with
missing information were not removed from the data set.

Figure 1: Averaged feature importances over 30 trials. Physics GRE score and
undergraduate GPA, appearing in orange, were the only factors found to be
meaningful.

##  III Results

Across the 30 runs, the average accuracy of our model predicting on the held-
out data was $75.8\%\pm 0.6\%$ and the average AUC was $.821\pm.002$. As our
model’s accuracy is significantly higher than the null accuracy of $52.3\%$,
the percent of students who were not accepted, and our AUC is above 0.7, our
model can be considered a reasonable model of the data. The accuracy and AUC
were also stable with regard to how many trees we used to construct the model
and the fraction of the data we used to train the model, with the accuracy
changing by at most 1% and the AUC by 0.005.

The feature importances averaged over the 30 runs are shown in Fig. 1. We find
numerical factors such as the applicant’s score on the physics GRE, the
applicant’s undergraduate GPA, the applicant’s quantitative and verbal GRE
scores, and their proposed research area to be more important in the
application process than any factor describing the applicant’s undergraduate
institution. Using recursive backward elimination to determine the meaningful
factors, we find the applicant’s physics GRE score and their undergraduate GPA
to be the only meaningful factors.

To verify that the applicant’s physics GRE score and their undergraduate GPA
were indeed the only meaningful factors, we then reran our random forest model
30 times using only these two factors as the predictors. Our average accuracy
was then $75.8\%\pm 0.6\%$ and our average area under the curve was
$.817\pm.003$, which are not statistically different from the values we found
using all fourteen factors shown in table 1. Further, if these factors are
truly meaningful, then plotting an applicant’s physics GRE score against their
undergraduate GPA should show distinct groups. To test this idea, we generated
all possible pairs of undergraduate GPA and physics GRE scores and ran them
through our model to find the predicted admissions decision. The results are
shown in Fig. 2 with the data from all of the applicants in our study
overlaid. While there does not appear to be a hard cutoff physics GRE score,
there appears to be a rough cutoff around 700. Interestingly, the decision
boundary has a negative slope between GPA of 3.6 and 4.0, suggesting a high
GPA can partially compensate for a physics GRE score below the soft threshold.

Figure 2: Plot of all applicant’s physics GRE scores vs their undergraduate
GPAs. The background coloring expresses the prediction of the model had an
applicant had that score.

##  IV Discussion

Perhaps unsurprisingly, we find numerical measures are the most important
factors for determining whether a domestic applicant will be accepted into
this physics graduate program, consistent with findings that graduate programs
with a large number of applicants use numerical measures as a first pass to
evaluate applicants [5, 3]. While we find no evidence of a minimum physics GRE
score, we do find evidence of a "rough cutoff" as described in Potvin et al.
around 700. Nevertheless, some students who scored significantly above this
threshold were not admitted. While we do not know the reasons why these
students were not admitted, Posselt noted that faculty may not admit superior
applicants if they do not believe the applicant will actually enroll in their
program [5].

Overall, our findings of the meaningful factors for admission to a physics
graduate program are consistent with Potvin et al.’s findings obtained by
surveying physics graduate admissions directors. Notably, we also find that
the physics GRE score, undergraduate GPA, quantitative GRE score, and proposed
research area are more important than other factors while the undergraduate
institution, GRE written score, and proximity/familiarity are less important
factors.

Interestingly, the program studied here appears to place more emphasis on the
verbal GRE score than the average program. This may be because our study only
looked at domestic students while Potvin et al.’s looked at all applicants.
Because international students also take the TOEFL while domestic students do
not and admissions directors ranked the TOEFL as more important the verbal
GRE, the TOEFL may take the place of the verbal GRE and hence lower the
perceived value of the verbal GRE relative to other factors.

Despite prior work suggesting institutional characteristics play an important
role in graduate admissions, we did not find institutional or departmental
characteristics to be meaningful to our model. Our result could be due to
differences in methodology or due to institutional effects being influential
but not dominant factors [22]. Indeed, Posselt suggests institutional factors
may be used to differentiate applicants with similar GPAs and GRE scores [5].
Therefore, we may not have found institutional factors to be meaningful since
they are used when primary factors such as GPA and physics GRE scores do not
sufficiently separate applicants.

While we did not have access to other criterion included in the Potvin et al.
such as application essays, research experiences, and recommendation letters,
we were still able to create a model that correctly predicted whether an
applicant would be admitted with 75% accuracy based solely off the applicant’s
undergraduate GPA and physics GRE score. While undergraduate GPA is a
significant predictor of completing a physics PhD, physics GRE is not, as
those scoring near the top of the physics GRE only have a 7% higher
probability of completing their PhD than those scoring near the bottom [27].
As the GRE is not associated with completing a doctoral degree and is known to
favor persons from majority groups in science [28], the outsized role of the
GRE in the admissions process should be questioned. Indeed, the American
Association of Physics Teachers recently released a recommendation against
using the physics GRE in graduate admissions [29].

##  V Limitations

As with any model, there are a few limitations to our model. First, the data
we used to make our model was not all the data that would be available to a
faculty member evaluating an application. Therefore, it is possible that
meaningful features other than GPA and physics GRE score could lie in the data
that was unavailable to us.

Second, feature importances can be inflated when features used in the model
are correlated, possibly causing some factors to be considered more important
than they actually are. There do exist methods for controlling for
correlations in random forest; however, unlike the conditional inference
forest algorithm used here, those methods can only be used with data
containing no missing values. Our initial analysis controlling for
correlations suggests that undergraduate GPA and physics GRE scores are still
the meaningful factors.

##  VI Future Work and Conclusion

Our work adds to the broader literature about graduate admissions and the
process by which applicants are judged. Because minority students may not
apply to graduate programs if they do not think they will be accepted,
elucidating the factors that determine whether an applicant will be accepted
is crucial. Simply increasing the number of applicants from underrepresented
groups will not increase their representation unless corresponding efforts are
made to admit these students. While our results align with prior work, these
results represent only one institution and may not be representative of all
United States physics PhD programs. Our future work will apply a similar
analysis to other physics graduate programs to test the generalizability of
our results. Given the unique structure of the admissions process at this
university, graduate programs with a more traditional admissions process may
assign different weights to the various parts of an application.

Further, the university studied here has recently moved to a rubric-based
admissions format, designed to take into account non-cognitive competencies
and program fit in addition to the more traditional admissions criteria such
as GPA and GRE scores. Our future work will examine how including these new
criteria may change the factors that are most predictive of an applicant being
admitted to the program.

###### Acknowledgements.

We would like to thank the graduate program directors at the studied
university for providing application data and for providing insight into their
admissions process.

## References

  * Ivie [2018] R. Ivie, Beyond Representation: Data to Improve the Situation of Women and Minorities in Physics and Astronomy (2018).
  * Cochran  _et al._ [2018] G. L. Cochran, T. Hodapp, and E. E. A. Brown, Identifying barriers to ethnic/racial minority students’ participation in graduate physics, in  _Physics Education Research Conference Proceedings_, PER Conference (Cincinnati, OH, 2018) pp. 92–95.
  * Potvin  _et al._ [2017] G. Potvin, D. Chari, and T. Hodapp, Investigating approaches to diversity in a national survey of physics doctoral degree programs: The graduate admissions landscape, Physical Review Physics Education Research 13, 020142 (2017).
  * Chari and Potvin [2019a] D. Chari and G. Potvin, Admissions practices in terminal master’s degree-granting physics departments: A comparative analysis, Physical Review Physics Education Research 15, 010104 (2019a).
  * Posselt [2016] J. R. Posselt,  _Inside graduate admissions_ (Harvard University Press, 2016).
  * Chari and Potvin [2019b] D. Chari and G. Potvin, Understanding the importance of graduate admissions criteria according to prospective graduate students, Physical Review Physics Education Research 15, 023101 (2019b).
  * [7] J. Posselt, T. Hernandez, G. Cochran, and C. Miller, Metrics First, Diversity Later? Making the Shortlist and Getting Admitted to Physics PhD Programs, Journal of Women and Minorities in Science and Engineering 10.1615/JWomenMinorScienEng.2019027863.
  * Gupta  _et al._ [2016] N. Gupta, A. Sawhney, and D. Roth, Will I Get in? Modeling the Graduate Admission Process for American Universities, in  _2016 IEEE 16th International Conference on Data Mining Workshops (ICDMW)_ (2016) pp. 631–638.
  * Waters and Miikkulainen [2014] A. Waters and R. Miikkulainen, GRADE: Machine Learning Support for Graduate Admissions, AI Magazine 35, 64 (2014).
  * Lux  _et al._ [2016] T. Lux, R. Pittman, M. Shende, and A. Shende, Applications of Supervised Learning Techniques on Undergraduate Admissions Data, in  _Proceedings of the ACM International Conference on Computing Frontiers_, CF ’16 (ACM, New York, NY, USA, 2016) pp. 412–417, event-place: Como, Italy.
  * Basu  _et al._ [2019] K. Basu, T. Basu, R. Buckmire, and N. Lal, Predictive Models of Student College Commitment Decisions Using Machine Learning, Data 4, 65 (2019).
  * Moore [1998] J. Moore, An expert system approach to graduate school admission decisions and academic performance prediction, Omega 26, 659 (1998).
  * Posselt [2015] J. R. Posselt, Disciplinary Logics in Doctoral Admissions: Understanding Patterns of Faculty Evaluation, The Journal of Higher Education 86, 807 (2015).
  * [14] I. U. C. for Postsecondary Research,  _The Carnegie Classification of Institutions of Higher Education, 2018 edition_ (Bloomington, IN).
  * Nicholson and Mulvey [2016] S. Nicholson and P. J. Mulvey,  _Roster of Physics Departments with Enrollment and Degree Data, 2015_, Tech. Rep. (American Institute of Physics, 2016).
  * Paxton and Bollen [2003] P. Paxton and K. A. Bollen, Perceived Quality and Methodology in Graduate Department Ratings: Sociology, Political Science, and Economics, Sociology of Education 76, 71 (2003).
  * Strobl  _et al._ [2007] C. Strobl, A.-L. Boulesteix, A. Zeileis, and T. Hothorn, Bias in random forest variable importance measures: Illustrations, sources and a solution, BMC Bioinformatics 8, 25 (2007).
  * Araújo  _et al._ [2005] M. B. Araújo, R. G. Pearson, W. Thuiller, and M. Erhard, Validation of species climate impact models under climate change, Global Change Biology 11, 1504 (2005).
  * Janitza  _et al._ [2013] S. Janitza, C. Strobl, and A.-L. Boulesteix, An AUC-based permutation variable importance measure for random forests, BMC Bioinformatics 14, 119 (2013).
  * Díaz-Uriarte and Alvarez de Andrés [2006] R. Díaz-Uriarte and S. Alvarez de Andrés, Gene selection and classification of microarray data using random forest, BMC Bioinformatics 7, 3 (2006).
  * Young  _et al._ [2019] N. T. Young, G. Allen, J. M. Aiken, R. Henderson, and M. D. Caballero, Identifying features predictive of faculty integrating computation into physics courses, Physical Review Physics Education Research 15, 010114 (2019).
  * Attiyeh and Attiyeh [1997] G. Attiyeh and R. Attiyeh, Testing for bias in graduate school admissions, The Journal of Human Resources; Madison 32, 524 (1997).
  * R Core Team [2018] R Core Team,  _R: A Language and Environment for Statistical Computing_ (R Foundation for Statistical Computing, Vienna, Austria, 2018).
  * Hothorn  _et al._ [2006] T. Hothorn, P. Bühlmann, S. Dudoit, A. Molinaro, and M. J. Van Der Laan, Survival ensembles, Biostatistics 7, 355 (2006).
  * Strobl  _et al._ [2008] C. Strobl, A.-L. Boulesteix, T. Kneib, T. Augustin, and A. Zeileis, Conditional variable importance for random forests, BMC Bioinformatics 9, 307 (2008).
  * Hapfelmeier  _et al._ [2014] A. Hapfelmeier, T. Hothorn, K. Ulm, and C. Strobl, A new variable importance measure for random forests with missing data, Statistics and Computing 24, 21 (2014).
  * Miller  _et al._ [2019] C. W. Miller, B. M. Zwickl, J. R. Posselt, R. T. Silvestrini, and T. Hodapp, Typical physics Ph.D. admissions criteria limit access to underrepresented groups but fail to predict doctoral completion, Science Advances 5, eaat7550 (2019).
  * Miller and Stassun [2014] C. Miller and K. Stassun, A test that fails, Nature 510, 303 (2014).
  * noa [2019] Statement on the Use of the GRE in Admissions to Graduate Physics Programs, American Association of Physics Teacher (2019).