11institutetext: Zentrum für Astronomie der Universität Heidelberg,
Landessternwarte, Königstuhl 12, 69117 Heidelberg, Germany
11email<EMAIL_ADDRESS>22institutetext: Department of
Astronomy, Stockholm University, AlbaNova University Center, SE-106 91
Stockholm, Sweden
33institutetext: Research School of Astronomy and Astrophysics, Australian
National University, Canberra, ACT 2611, Australia
44institutetext: ARC Centre of Excellence for All Sky Astrophysics in 3
Dimensions (ASTRO 3D), Australia

# A method for identifying metal-poor stars with
Gaia BP/RP spectra ††thanks: Table A.1 and a table with all the coefficients
of the method are only available in electronic form at the CDS via anonymous
ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via
http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/

T. Xylakis-Dornbusch 11Member of the International Max Planck Research School
for Astronomy & Cosmic Physics at the University of Heidelberg (IMPRS-
HD).Member of the International Max Planck Research School for Astronomy &
Cosmic Physics at the University of Heidelberg (IMPRS-HD). N. Christlieb 11
K.Lind 22 T.Nordlander 3344

(Received )

###### Abstract

Context. The study of the oldest and most metal-poor stars in our Galaxy
promotes our understanding of the Galactic chemical evolution and the
beginning of Galaxy and star formation. However, they are notoriously
difficult to find, with only five stars at $\mathrm{[Fe/H]<-5.0}$ having been
detected to date. Thus, the spectrophotometric data of 219 million sources
which became available in the third Gaia Data Release comprise a very
promising dataset for the identification of metal-poor stars.

Aims. We want to use the low-resolution Gaia Blue Photometer / Red Photometer
(BP/RP) spectra to identify metal-poor stars. Our primary aspiration is to
help populate the poorly constrained tail of the metallicity distribution
function of the stellar halo of the Galaxy.

Methods. We developed a metal-poor candidate selection method based on flux
ratios from the BP/RP Gaia spectra, using simulated synthetic spectra.

Results. We found a relation between the relative iron abundance and the flux
ratio of the Ca H & K region to that of the $\mathrm{H\beta}$ line. This
relation is temperature and surface gravity dependent, and it holds for stars
with $\mathrm{4800\,K\leq T_{eff}\leq 6300\,K}$. We applied it to noisy
simulated synthetic spectra and inferred $\mathrm{[Fe/H]}$ with an uncertainty
of $\sigma_{\mathrm{[Fe/H]}}\lessapprox 0.65$ dex for
$\mathrm{-3\leq[Fe/H]}\leq 0.5$ and G=15-17mag, which is sufficient to
identify stars at $\mathrm{[Fe/H]<-2.0}$ reliably. We predict that by
selecting stars with inferred $\mathrm{[Fe/H]}\leq-2.5$ dex, we can retrieve
80% of the stars with $\mathrm{[Fe/H]}\leq-3$ and have a success rate of about
50%, that is one in two stars we select would have $\mathrm{[Fe/H]}\leq-3$. We
do not take into account the effect of reddening, so our method should only be
applied to stars which are located in regions of low extinction.

###### Key Words.:

stars: Population II – surveys – stars: carbon

##  1 Introduction

The Gaia survey (Gaia Collaboration et al., 2016) poses an unprecedented
opportunity to understand the structure, history, and evolution of our Galaxy
better. Aside from the astrometry and photometry for over a billion stars,
Gaia has also delivered spectrophotometric data in its third Data Release
(Gaia DR3), which was made public on June 13, 2022. This could serve as a tool
to map out the Galaxy in terms of the relative iron abundance
$\mathrm{[Fe/H]}$ – which is a proxy for the total amount of metals in a star
– and consequently find new metal-poor stars. Stars deficient in elements
heavier than helium, called metal-poor stars, are usually very old, with ages
comparable to that of the Universe (Frebel & Norris, 2013). The basic idea is
that stellar atmospheres preserve, to a large extent, the makeup of their
birth cloud, hence a metal deficient star should have formed along with the
first few generations of stars after the Big Bang. That is why these objects
are very interesting, and they can facilitate our understanding of the past,
even though they can be found close by. In particular, they can be used as
tracers of the evolution of the elements and of the Galaxy, help us understand
how the distribution of the first population of stars (Pop III) looked, and
assist in constraining Big Bang nucleosynthesis further (Frebel, 2010; Beers &
Christlieb, 2005; Sneden et al., 2008).
The Gaia Collaboration itself, through the Data and Analysis Consortium
(DPAC), has already, and will also in the future, deliver astrophysical
parameters using the Gaia data. With the astrometric and photometric data of
Gaia DR2, Andrae et al. (2018) delivered temperatures, extinctions and
reddening. DPAC also used the Gaia DR3 Blue Photometer / Red Photometer
(BP/RP) spectra to estimate metallicities of stellar objects (Creevey et al.,
2022; Fouesneau et al., 2022; Andrae et al., 2022). Liu et al. (2012)
predicted that Aeneas – which is a Bayesian method employed by DPAC for the
stellar parameter estimation – would estimate the metallicity with an accuracy
of 0.1-0.2 dex, given that the extinction is $\mathrm{A_{0}<2\,mag}$ and that
the magnitude is $G=15$mag. Now with Gaia DR3 available, Andrae et al. (2022)
estimated $\mathrm{[M/H]}$ values for several million sources, and provide
information on how to use them.
The search for metal-poor stars is as pertinent as ever now. The ongoing
photometric SkyMapper Southern Sky (SMSS) survey is actively searching for
extremely metal-poor stars ($\mathrm{[Fe/H]<-3\,dex}$ according to the
classification of Beers & Christlieb (2005)). With its DR1.1, the SkyMapper
team found that 40% of their candidates had $\mathrm{[Fe/H]<-2.75\,dex}$ and
only $\approx 7\%$ had $\mathrm{[Fe/H]>-2\,dex}$, that is an efficacy of
$\approx 93\%$ in finding stars with $\mathrm{[Fe/H]<-2\,dex}$ (Da Costa et
al., 2019). SkyMapper also led to the discovery of the star with the lowest
detected Fe abundance (Nordlander et al., 2019). The Pristine survey, which is
being carried out in the northern hemisphere, was tailored to find metal-poor
stars with a narrow-band filter centered on the metallicity-sensitive Ca II H
& K lines (Starkenburg et al., 2017). Pristine has a 56% reported efficacy of
finding stars with $\mathrm{[Fe/H]<-2.5\,dex}$ and 23% for stars with
$\mathrm{[Fe/H]<-3\,dex}$ (Aguado et al., 2019). The objective of this work is
to develop a recipe that can efficiently identify metal-poor stars, and
especially stars within and below the extremely metal-poor regime
($\mathrm{[Fe/H]<-3\,dex}$). In Section 2 we present the reasoning of our
approach for this endeavor, as well as the tools we used. We also describe the
parameter space where our method is applicable, along with a detailed
description of the development of our technique. Further, in Section 3 we
present our results from applying our method on noise-free and noisy spectra.
Therein, we also investigate the dependence of our procedure on the different
astrophysical parameters, and we additionally study the effect of extinction.
Lastly, we explore the expected efficiency of our technique.

##  2 Methods

For this work, we used the Ulysses Simulator (Astraatmadja, 2015) and
synthetic spectra (Nordlander et al., 2019) (see Sections 2.1.2 and 2.1.1,
respectively) in order to simulate the spectrophotometric data of BP and RP,
respectively, on board Gaia. The very low resolution of the expected BP/RP
spectra (see Figure 1) inclined us to use integrated fluxes of different parts
of the spectra, such that a relation varying with metallicity could be found.
Those spectral areas had to be sensitive to the change in metallicity and also
to at least one other quantity that can be known a priori or that can be
directly derived from the spectra themselves.

Figure 1: Synthetic spectrum before (orange) and after (green and magenta)
simulation. The spectrum has $\mathrm{T_{eff}=5300\,K,\,\log
g=2.5\,dex,\,[Fe/H]=-2\,dex}$ and $\mathrm{[C/Fe]=0.0\,dex}$.

###  2.1 Data and simulations

Choosing the temperature and surface gravity space for our dataset was the
first critical step for our method. The desired relation has to hold in that
parameter space, which would then allow us to use those parameters as priors
when applying the method described in this paper. Additionally, the dataset
has to cover a wide metallicity range, so that both metal-poor and metal-rich
stars can be covered, and we can ultimately distinguish between them. Since
low metallicities are involved in the process,, the last important factor that
we have to pay attention to is the carbon enhancement that has been observed
in many metal-poor stars (Beers & Christlieb, 2005; Frebel & Norris, 2015;
Arentsen et al., 2021). The initial dataset parameters are shown in Table 1
and Figure 2.

Table 1: Initial parameter space of the synthetic spectra.

Parameter | Range | Step
---|---|---
$\mathrm{T_{eff}}$ (K) | $3500\,\leq\mathrm{T_{eff}}\leq 7800$ | 100
$\mathrm{\log g}$ (dex) | $0.5\leq\mathrm{\log g}\leq 5.0$ | 0.5
$\mathrm{[Fe/H]}$ (dex) | $-7.0\leq\mathrm{[Fe/H]}\leq 0.5$ | 0.5
$\mathrm{[C/Fe]}$ (dex) | $-1.5\leq\mathrm{[C/Fe]}\leq 6.0$ | 0.5
111The metallicity and carbon enhancement combinations are restricted to
physically meaningful combinations (see Figure 2).

Figure 2: Parameter space we used to develop our method. Top plot: Each point
in the top plot represents a specific combination of relative carbon and iron
abundances. Our entire dataset is comprised of all the meaningful temperature
and surface gravity combinations (see bottom plot) for each and every one of
these relative abundance combinations. Bottom plot: A PARSEC isochrone of 12
Gyr and $\mathrm{[Fe/H]=-3.0\,dex}$, which we used to determine the
astrophysical parameters of our preliminary dataset. Overplotted are the
temperature-surface gravity pairs of the final parameter space we used to
develop our method.

####  2.1.1 Synthetic spectra

For this work, we used synthetic spectra from a one dimensional (1D) grid
where local thermodynamic equilibrium (LTE) is assumed (Nordlander et al.,
2019), that was computed with the Turbospectrum code (Alvarez & Plez, 1998;
Plez, 2012) and MARCS model atmospheres (Gustafsson et al., 2008). Nordlander
et al. (2019) adopted the Solar chemical composition from Asplund et al.
(2009); for $\mathrm{[Fe/H]<-1.0}$, $\mathrm{[\alpha/Fe]}=+0.4$ was adopted;
for models with $\mathrm{-1.0<[Fe/H]<0.0}$, a relationship of
$\mathrm{[\alpha/Fe]=+0.4\cdot[Fe/H]}$ was used. The models take continuum
scattering into account and use $v_{mic}=1\,\mathrm{kms^{-1}}$ for models with
$\mathrm{\log g}\geq 4.0$. Models with $\mathrm{\log g}\leq 3.5$ take
spherical symmetry into account and assume $v_{mic}=2\,\mathrm{kms^{-1}}$. The
line data are a combination of atomic data from VALD3 (Ryabchikova et al.,
2015) and molecular data from a variety of sources; here the primary molecular
line list is the one for CH from Masseron et al. (2014). This grid also
contains spectra with varying carbon abundances, which we used for this work.

####  2.1.2 Ulysses

Ulysses (Astraatmadja, 2015) is a simulator which takes spectra as input and
delivers the equivalent BP/RP spectrum as it would be observed from Gaia. The
final DR3 BP/RP spectra will be a combination of many different epoch
observations, spanning up to the entire time of the missions’ operation. That
is why the input spectrum is being convolved with an averaged line spread
function (LSF). For each of the two spectrophotometers (BP and RP), an
averaged LSF exists, which is the product of all the LSFs of the different
CCDs and telescopes. The parameters of Ulysses that can be tuned, among
others, are the number of transits – which are the number of times an object
was observed – the amount of interstellar extinction $A_{0}$, the apparent
magnitude of the source, as well as the spectrum oversampling. The spectrum
oversampling is essentially the number of spectra that will be combined into
the final BP/RP spectrum. Each spectrum is sampled over 60 pixels, but for
different observations, different parts of the spectrum are being sampled, so
that there will be more than 60 pixels of data in the end. Carrasco et al.
(2021) present a model for the internal calibration of the BP/RP spectra, in
particular how to produce a mean spectrum from all the epoch spectra of the
same source. The extinction curve models that Ulysses implements are those of
Cardelli et al. (1989) and Fitzpatrick (1999). We selected the Fitzpatrick
(1999) models for our simulations, after determining both models produce very
similar results.
The core products of Ulysses are noise-free BP/RP spectra, noise-free Gaia
photometry, and extinctions. Other outputs are also possible, for instance the
end-of-mission noisy spectra, which are being generated with the addition of
noise to the noise-free spectra. The noise model employed by Ulysses is the
one from Jordi et al. (2010). An overview of the configuration we used to
simulate our data is shown in Table 2.

####  2.1.3 Parameter space

A $\chi^{2}$ test of the simulated spectra showed that they carry enough
information in the region of the Ca II H & K lines to distinguish between
metal-rich stars and stars with $\mathrm{[Fe/H]<-2.0}$. However, this was
observed only for stars with $\mathrm{T_{eff}\geq 4800\,K}$ and also depending
on their surface gravity. A $\mathrm{12\,Gyr}$ and $\mathrm{Z=0.001}$
isochrone (Bressan et al., 2012) was used to choose realistic
$\mathrm{T_{eff}}$ \- $\mathrm{\log g}$ pairs for the preliminary dataset (see
Figure 2).

Table 2: Configuration of the Ulysses simulator which we used in order to
construct our method.

Parameter | Value
---|---
G magnitude (mag) | 15
spectrum oversampling | 4
number of transits | 75
extinction $\mathrm{A_{0}}$ | 0.0
222The method is described in Section 2.1.2. The same configuration was used
for the test and applications described in Section 3, apart from the magnitude
and extinction which are explicitly mentioned when they have different values
than in the table above.

###  2.2 Flux ratios

Differential simulated spectra were used to find out how the flux changes with
changing metallicity. For this exercise, we used the spectra with
$\mathrm{[Fe/H]=-2.5\,dex}$ as a reference, and we subtracted the spectra of
lower and higher metallicity that had all the other astrophysical parameters
identical to the reference spectrum (Figure 3). The regions we expected to be
$\mathrm{[Fe/H]}$ sensitive were the Ca II H & K and G-band wavelength ranges,
which was confirmed: the lower the metallicity of the object, the higher the
flux. The regions of the spectra that cover the Ca II near-infrared triplet as
well as the H$\beta$ absorption line showed a greater sensitivity to the
change in all the other stellar parameters compared to $\mathrm{[Fe/H]}$.
Using the ratios of the aforementioned spectral regions, that is the ratio of
the integrated Ca H & K flux to that of the H$\beta$ region
($fr\mathrm{{}_{CaHK/H\beta}}$), and that of the G band integrated flux to the
Ca NIR triplet ($fr\mathrm{{}_{G/CaNIR}}$), we can see a relation with
metallicity (Figure 4). Figure 4 was created from our entire dataset, that is
temperature, surface gravity (Figure 2), and all physically meaningful
$\mathrm{[C/Fe]}$-$\mathrm{[Fe/H]}$ combinations, in other words
$\mathrm{[C/Fe]+[Fe/H]\leq-1\,dex}$ for $\mathrm{[Fe/H]}\leq 0.5$ and
$-1.5\geq\mathrm{[C/Fe]\leq 0}$ for $\mathrm{[Fe/H]}\geq-0.5$.

Figure 3: Differential BP (top panel) and RP (bottom panel) spectra of Solar
$\mathrm{[C/Fe]}$-scaled stars, with different $\mathrm{[Fe/H]}$. Each
differential spectrum results from $reference\,spectrum-spectrum,$ where
$\mathrm{[Fe/H]_{ref}=-2.5}$. The shaded areas correspond to the regions we
used for our flux ratios, i.e., the Ca II H & K region (purple), the G-band
(magenta), the H$\beta$ (yellow), and the Ca II NIR (green) regions.

###  2.3 Metallicity-temperature-$fr\mathrm{{}_{CaHK/H\beta}}$ relation

The relation between $fr\mathrm{{}_{CaHK/H\beta}}$, $fr\mathrm{{}_{G/CaNIR}}$,
and $\mathrm{[Fe/H]}$ that is emerging in Figure 4 encouraged us to try to
find a relation of the form
$F(fr\mathrm{{}_{CaHK/H\beta}},fr\mathrm{{}_{G/CaNIR}})=\mathrm{[Fe/H]}$. In
this venture, we found that $fr\mathrm{{}_{G/CaNIR}}$ is – by itself – an
indicator for the effective temperature (see Figure 5), but it does not carry
any information about the metallicity. On the other hand, for constant
$fr\mathrm{{}_{G/CaNIR}}$, the metallicity is changing smoothly with the
change of $fr\mathrm{{}_{CaHK/H\beta}}$ (Figure 6), hence
$fr\mathrm{{}_{CaHK/H\beta}}$ carries information concerning the iron
abundance. There is, however, a scatter that accounts for variations in
$\mathrm{T_{eff}}\,,\mathrm{\log g}$, and $\mathrm{[C/Fe]}$. So, since
$fr\mathrm{{}_{G/CaNIR}}$ is not metallicity-sensitive, as mentioned above,
but rather a temperature indicator, we thought it best to proceed with the
objective of developing a relation which is comprised of
$fr\mathrm{{}_{CaHK/H\beta}}$ and some relevant stellar parameters. The usage
of stellar parameters makes the range of applicability more straightforward to
implement, and subsequently more user-friendly. By inspecting Figure 6, we
assessed that the therein $fr\mathrm{{}_{CaHK/H\beta}}$-$\mathrm{[Fe/H]}$
correlation can be described with an exponential relation (Figures 4 and 6),
which was confirmed by the residuals of the fit. This exponential relation can
then be expressed with respect to the temperature and surface gravity of the
model stars (Figure 12), enabling us to use them as priors, for example from
Gaia itself (Andrae et al., 2018). Finally, we can use the
iso-$fr\mathrm{{}_{CaHK/H\beta}}$ lines (see Figure 12), as well as the
effective temperatures and surface gravities, in order to infer the
metallicity:

$F(fr\mathrm{{}_{CaHK/H\beta}},\mathrm{T_{eff}},\mathrm{\log
g})=\mathrm{[Fe/H]}$ (1)

$\displaystyle F=-(\mathrm{T_{eff}}\cdot\mathrm{\log g})\cdot e^{b\cdot
fr\mathrm{{}_{CaHK/H\beta}}+c}+d,$ (2)

where $b,c,$ and $d$ are $\mathrm{T_{eff}}$ and $\mathrm{\log g}$ dependent
coefficients.

Figure 4: Smoothly changing metallicity for model stars of G=15 mag. On the
inset plot, it is shown how well the noisy simulated flux ratios track the
modeled (noiseless) ones. Figure 5: $fr\mathrm{{}_{G/CaNIR}}$ is an effective
temperature indicator; however, it does not carry information regarding the
iron abundance. Figure 6: Exponential decline of $fr\mathrm{{}_{CaHK/H\beta}}$
with increasing [Fe/H], for a roughly constant $fr\mathrm{{}_{G/CaNIR}}$. This
behavior starts to break down for $\mathrm{T_{eff}<4800\,K}$. The scatter at
each iron abundance reflects the different carbon enhancement and surface
gravity, as well as a moderate ($\sim 300$K) variation in temperature.

The coefficients $b,c,d$ are the result of fitting Eq. 2 to
$fr\mathrm{{}_{CaHK/H\beta}}$-metallicity pairs of a roughly constant
temperature – $\mathrm{T_{eff}}\mp 30$K – and constant surface gravity (see
Figure 13) for metallicities greater than $-3.5$ dex. The aforementioned
$fr\mathrm{{}_{CaHK/H\beta}}$-$\mathrm{[Fe/H]}$ pairs result from the
iso-$fr\mathrm{{}_{CaHK/H\beta}}$ lines.

##  3 Results

We first tested our method on the very same data we used to construct it, and
then we examined how it is influenced by different stellar parameters.
Finally, we applied our method to noisy spectra.

###  3.1 Testing on the model spectra

We applied our method on the same model spectra we used to find relation 2.
Figure 7 shows that our method works very well for stars with
$\mathrm{[Fe/H]}\geq-3$ dex since $\sigma_{\mathrm{{[Fe/H]}}}\leq 0.6$ dex.
Specifically, for red giant branch stars $\sigma_{\mathrm{{[Fe/H]}}}=0.6$ and
for turnoff stars, $\sigma_{\mathrm{{[Fe/H]}}}\approx 0.3$ at
$\mathrm{[Fe/H]}=-3$. Further, even for stars with $\mathrm{[Fe/H]}=-3.5$ dex,
$\sigma$ is smaller than $1.0$ dex. On the other hand, we can identify stars
with $\mathrm{[Fe/H]}\geq-1$ very accurately, that is to say with
$\sigma_{\mathrm{{[Fe/H]}}}\approx 0.1$ dex. The bias of these results (see
Figure 7), which we defined as
$\overline{\mathrm{[Fe/H]_{inf}}-\mathrm{[Fe/H]_{ref}}}$, shows that our
method tends to overestimate the metallicity for
$\mathrm{[Fe/H]_{ref}}\leq-1$, in other words it assigns higher values than
the true ones. The benefit of this result is that the inferred metallicities
that are in and below the metal-poor range are most probably as low as or even
lower than the true ones (see Section 3.3). Further, the way the difference
between the bias and $\sigma$ decreases as we go to lower
$\mathrm{[Fe/H]_{ref}}$ indicates that there is a point below which we cannot
distinguish metallicities, and that is around $\mathrm{[Fe/H]_{ref}}=-3.5$
dex.

###  3.2 Application on noisy spectra and the dependence on stellar parameters

We subsequently applied our method to noisy spectra. For each of our model
spectra (G=15mag), we generated 20 noisy ones: we inferred the metallicity
with $\sigma_{\mathrm{{[Fe/H]}}}\leq 0.6$ dex for $\mathrm{{[Fe/H]}}\geq-3$
dex (Figure 8). Again, the uncertainty in inferring metallicity for metal-poor
stars ($\mathrm{[Fe/H]}=-1$) and above is very low with
$\sigma_{\mathrm{{[Fe/H]}}}\lessapprox 0.12$ dex. Further we sought to
determine how our method performs in the different temperature bins. For this
purpose, we computed $\sigma_{\mathrm{{[Fe/H]}}}$ for all spectra of the same
temperature, surface gravity, and relative carbon abundance (Figure 14). We
found that the dependence on temperature is linked with $\mathrm{[C/Fe]}$.
When $\mathrm{[C/Fe]}\leq 0.0$ dex, our method performs slightly better for
lower effective temperatures of the same $\mathrm{\log g}$ when
$\mathrm{{[Fe/H]}<-4}$. In contrast, when the iron abundance is greater than
-4, performance becomes independent of $\mathrm{T_{eff}}$. When
$\mathrm{[C/Fe]>0}$, performance has a similar behavior for
$\mathrm{{[Fe/H]}<-4}$, that is to say it is better for lower temperatures.
However, above that threshold, the effect is reversed: the performance is much
better for higher effective temperatures. That phenomenon is intensified as
$\mathrm{[C/Fe]}$ grows. We expected that we could determine the metallicity
for lower temperatures more accurately since the Ca H & K lines are stronger
in cooler stars compared to hotter ones, even down to very low metallicities.
This, as stated above, is the case for $\mathrm{[C/Fe]\leq 0}$, but not above
that threshold. The reason for this is probably due to the many carbon lines
that arise as a result of the high carbon enhancements, which in turn are also
stronger at lower effective temperatures and, consequently, our method assigns
higher metallicities to those spectra. Furthermore, we investigated the
dependence on different relative carbon to iron abundances, for which we
already have some information as stated above. To test that, we separated our
data into surface gravity and $\mathrm{[C/Fe]}$ bins, and calculated the error
of the inferred metallicities (Figure 15). There is an obvious difference in
the performance pattern for stars below and above $\mathrm{[C/Fe]}=0$ dex. At
and below that threshold, $\sigma_{\mathrm{[Fe/H]}}$ is almost independent of
$\mathrm{[C/Fe]}$. Above it, however, the performance declines as the relative
carbon abundance is rising. This effect though is attenuated as temperature
rises.
Another factor we considered is the extinction $A_{0}$. We used synthetic
spectra with $\mathrm{T_{eff}=5500}$ K, $\mathrm{\log g=3.5}$ dex,
$\mathrm{[Fe/H]=-2}$ dex, and all possible $\mathrm{[C/Fe]}$ combinations –
for which $\mathrm{[Fe/H]+[C/Fe]}\geq-1$ is true – and passed them through
Ulysses generating spectra of two kinds for each model: a noiseless spectrum
and 20 noisy ones. We expected that the performance would decline with greater
extinction, which is what we observed (see Figure 9). We validated this result
by performing the above exercise for a collection of spectra with varying
astrophysical parameters. As extinction rises, $fr\mathrm{{}_{CaHK/H\beta}}$
and $fr\mathrm{{}_{G/CaNIR}}$ decrease, in other words the model star appears
to be colder and the Ca H & K features are less distinguishable. When we use
Figure 4 as a reference, a model star with increasing extinction shifts left
and downwards in the plot. However, according to Christlieb et al. (2008),
about $90\%$ of stars at high galactic latitudes (i.e., at $|b|\geq 20^{o}$)
have a reddening that is smaller than 0.06 ($E(B-V)\leq 0.06$), which
corresponds to $A_{v}<0.19$. Thus, the issue of extinction is outside the
scope of this paper, but we consider addressing its effect in our follow-up
work.
Lastly, we wanted to find out for which magnitude our method starts to break
down. Hence, we generated noisy spectra (one noisy spectrum for each set of
stellar parameter combinations) of magnitudes G=16, 17, 18 mag. Up to 17 mag,
the performance of our method is roughly on the same levels of precision
(Figure 8), and at 18 mag, $\sigma_{\mathrm{[Fe/H]}}]<0.85$ for
$\mathrm{[Fe/H]\geq-3}$. Our method does not seem to break down, but it rather
gets less precise as the magnitude rises. Even so, at 18 mag, we can infer
metallicities down to $-2$ dex with an uncertainty of
$\sigma_{\mathrm{[Fe/H]}}]\leq 0.6$ dex.

Figure 7: $\sigma$ and bias of the inferred metallicities of our noise-free
spectra of G=15 mag. Figure 8: $\sigma$ and bias of the inferred metallicities
of noisy spectra of different magnitudes: G=15, 16, 17, 18 mag.

Figure 9: Testing the dependence of metallicity inference on different
extinction coefficients has the following expected outcome: the higher the
extinction, the higher the uncertainty. We used one set of stellar parameters
and metallicity, i.e., $\mathrm{T_{eff}=5500}$ K, $\mathrm{\log g=3.5}$ dex,
and $\mathrm{[Fe/H]=-2}$ dex, but various relative carbon abundances. On the
top are the $\sigma_{\mathrm{{[Fe/H]}}}$ and bias in bins of extinctions, and
in the bottom the results are also in bins of different $\mathrm{[C/Fe]}$.

###  3.3 Predicted success rate

Lastly, an application of our method to a simulated dataset (G=15 mag) of
realistic $\mathrm{[Fe/H]}$ and $\mathrm{[C/Fe]}$ distributions (Figure 10)
served the purpose of assessing the expected success rate of our method.
Additionally, this test enabled us to set the metal-poor threshold, that is to
say the derived $\mathrm{[Fe/H]}$ value below which metal-poor stars can be
selected, with the greatest efficiency and completeness. We constructed the
aforementioned dataset in the following manner: a) we randomly drew
$\mathrm{T_{eff}}$-$\mathrm{\log g}$-$\mathrm{[Fe/H]}$-$\mathrm{[C/Fe]}$
combinations from our parameter space; b) then we replaced the metallicity and
carbonicity values by drawing new ones from two different metallicity
distribution functions (MDF) of halo stars, and from two carbonicity
distribution functions, respectively. Specifically, we used the MDF from
Youakim et al. (2020) when the initially drawn metallicity was $>-2$ dex, and
the MDF from Placco et al. (2014) otherwise. Concerning the carbonicity, we
used the respective carbon-enhancement distribution for the Placco et al.
(2014) MDF (Placco et al., 2014), and we used the $\mathrm{[C/Fe]}$
distribution from Yoon et al. (2016) for $\mathrm{[Fe/H]}>-2$.
We found that when we selected all stars with an inferred
$\mathrm{[Fe/H]_{inf}}\leq-2.5$, we recovered $80\%$ of stars with
$\mathrm{[Fe/H]}\leq-3$, and we had a contamination of about 2% of stars with
metallicities above -2.5. It should be noted, however, that these
“contamination” stars all have $\mathrm{[Fe/H]}<-2$. Furthermore, about 55% of
the stars with $\mathrm{[Fe/H]_{inf}}\leq-2.5$ have a reference
$\mathrm{[Fe/H]}\leq-3$, which means that one in two of the selected stars
would at least be extremely metal-poor ($\mathrm{[Fe/H]}\leq-3$). If we were
to select the metallicity threshold at $\mathrm{[Fe/H]_{inf}}\leq-3$, the hit
rate for stars below -3 would increase to 9 out of 10. We would, however, fail
to detect about 75% of stars with $\mathrm{[Fe/H]}\leq-3$ (see Figure 11). The
overall results of this test are detailed in Figure 11. Specifically, we show
how the success rate, the completeness, and the contamination change by
selecting a different threshold, that is for
$\mathrm{[Fe/H]_{inf\,threshold}}=-2.0,-2.5,-3.0,-3.5$ dex. We define the
success rate as the percent of the selected stars that actually have
$\mathrm{[Fe/H]_{ref}}\leq-3.0$, the completeness as the percent of the total
number of stars with $\mathrm{[Fe/H]_{ref}}\leq-3.0$ that have
$\mathrm{[Fe/H]_{inf}}\leq\mathrm{[Fe/H]_{inf\,threshold}}$, and the
contamination as the percent of selected stars that have
$\mathrm{[Fe/H]_{ref}}>\mathrm{[Fe/H]_{inf}}$.

Figure 10: Metallicity and carbonicity of the simulated dataset we used to
predict the success rate of our method. The color bar designates the number of
models in each $\mathrm{[Fe/H]}$-$\mathrm{[C/Fe]}$ bin. Figure 11: Calculation
of the success rate, the completeness (left plot), and the contamination
(right plot) for threshold values $\mathrm{[Fe/H]_{inf}}=-3.5,-3.0,-2.5,-2.0$.
The contamination was calculated with reference to the threshold metallicity,
whereas the completeness and the success rate were calculated with reference
to model spectra with $\mathrm{[Fe/H]}\leq-3$.

##  4 Conclusions

We developed a method using flux ratios of metallicity-sensitive regions from
the Gaia BP/RP low resolution spectra in order to find new metal-poor stars.
This method is applicable when stars have $\mathrm{4800\,K\geq T_{eff}\leq
6300\,K}$. We took into account the fact that a large fraction of metal-poor
stars are carbon enhanced, and thus used a grid of synthetic spectra with
varying $\mathrm{[Fe/H]}$ and $\mathrm{[C/Fe]}$. We found an exponential
relation between the metallicity and the $fr\mathrm{{}_{CaHK/H\beta}}$ flux
ratio, which is temperature and surface gravity dependent. Therefore, our
method requires both of those stellar parameters as priors. We first applied
our method to the very same noiseless data ($\mathrm{G=15\,mag}$) we used to
construct it, and $\mathrm{[Fe/H]}$ was inferred with an uncertainty of
$\sigma\lessapprox 0.6$ dex at $\mathrm{[Fe/H]}\leq-3.0$ dex. Our method’s
performance was approximately on the same level for noisy spectra of the same
and greater magnitudes, that is $\sigma\lessapprox 0.65$ for
$\mathrm{[Fe/H]}\leq-3.0$ dex and $\mathrm{G=15,16,17\,mag}$. We found that
the performance depends on temperature, but in conjunction with the relative
carbon abundance: for Solar $\mathrm{[C/Fe]}$ and below, performance is
slightly better for lower temperatures of the same surface gravity when
determining $\mathrm{[Fe/H]}\lessapprox-3.5$. For $\mathrm{[Fe/H]}\geq-3$, the
performance is practically independent of $\mathrm{T_{eff}}$. When
$\mathrm{[C/Fe]}>0$, $\sigma_{\mathrm{[Fe/H]}\lessapprox-3.5}$ is as in the
Solar case, that is it is lower for lower $\mathrm{T_{eff}}$. However, the
determination of the metallicity above $-3.5$ dex presents a lower uncertainty
at higher temperatures of the same $\mathrm{\log g}$. A dependence of the
performance from the relative carbon-to-iron abundance is observed when
$\mathrm{[C/Fe]}>0$. Further, even for spectra of $\mathrm{18^{th}}$
magnitude, we can determine metallicities down to $-2$ with an uncertainty of
$\sigma_{\mathrm{[Fe/H]}}<0.6$ dex, which is sufficient for identifying metal-
poor stars. Nevertheless, we find that the overall performance of our method
deteriorates with rising extinction $A_{0}$. Lastly, we predict that by
selecting stars with $\mathrm{[Fe/H]_{inf}}\leq-2.5$ dex, we achieve a
completeness of 80% of the stars with $\mathrm{[Fe/H]}\leq-3$, and that one in
two stars we select is extremely metal-poor. Our method will be applied to
Gaia DR3 data and the results will be published in a forthcoming paper.

###### Acknowledgements.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”,
subproject A04). KL acknowledges funds from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(Grant agreement No. 852977) and funds from the Knut and Alice Wallenberg
foundation. This research was supported by the Australian Research Council
Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
through project number CE170100013. This work was supported by computational
resources provided by the Australian Government through the National
Computational Infrastructure (NCI) under the National Computational Merit
Allocation Scheme and the ANU Merit Allocation Scheme (project y89).

## References

  * Aguado et al. (2019) Aguado, D. S., Youakim, K., González Hernández, J. I., et al. 2019, Monthly Notices of the Royal Astronomical Society, 490, 2241–2253
  * Alvarez & Plez (1998) Alvarez, R. & Plez, B. 1998, A&A, 330, 1109
  * Andrae et al. (2018) Andrae, R., Fouesneau, M., Creevey, O., et al. 2018, A&A, 616, A8
  * Andrae et al. (2022) Andrae, R., Fouesneau, M., Sordo, R., et al. 2022, arXiv e-prints, arXiv:2206.06138
  * Arentsen et al. (2021) Arentsen, A., Starkenburg, E., Aguado, D. S., et al. 2021, Monthly Notices of the Royal Astronomical Society, 505, 1239–1253
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Astraatmadja (2015) Astraatmadja, T. L. 2015, Technical Report GAIA-C8-TN-MPIA-TLA-001, Max Planck Institute for Astronomy, Heidelberg
  * Beers & Christlieb (2005) Beers, T. C. & Christlieb, N. 2005, Annual Review of Astronomy and Astrophysics, 43, 531
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Carrasco et al. (2021) Carrasco, J. M., Weiler, M., Jordi, C., et al. 2021, A&A, 652, A86
  * Christlieb et al. (2008) Christlieb, N., Schörck, T., Frebel, A., et al. 2008, A&A, 484, 721
  * Creevey et al. (2022) Creevey, O. L., Sordo, R., Pailler, F., et al. 2022, arXiv e-prints, arXiv:2206.05864
  * Da Costa et al. (2019) Da Costa, G. S., Bessell, M. S., Mackey, A. D., et al. 2019, Monthly Notices of the Royal Astronomical Society, 489, 5900–5918
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63
  * Fouesneau et al. (2022) Fouesneau, M., Frémat, Y., Andrae, R., et al. 2022, arXiv e-prints, arXiv:2206.05992
  * Frebel (2010) Frebel, A. 2010, Astronomische Nachrichten, 331, 474
  * Frebel & Norris (2013) Frebel, A. & Norris, J. E. 2013, in Planets, Stars and Stellar Systems (Springer Netherlands), 55–114
  * Frebel & Norris (2015) Frebel, A. & Norris, J. E. 2015, Annual Review of Astronomy and Astrophysics, 53, 631–688
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951
  * Jordi et al. (2010) Jordi, C., Gebran, M., Carrasco, J. M., et al. 2010, A&A, 523, A48
  * Liu et al. (2012) Liu, C., Bailer-Jones, C. A. L., Sordo, R., et al. 2012, MNRAS, 426, 2463
  * Masseron et al. (2014) Masseron, T., Plez, B., Van Eck, S., et al. 2014, A&A, 571, A47
  * Nordlander et al. (2019) Nordlander, T., Bessell, M. S., Da Costa, G. S., et al. 2019, MNRAS, 488, L109
  * Placco et al. (2014) Placco, V. M., Frebel, A., Beers, T. C., & Stancliffe, R. J. 2014, ApJ, 797, 21
  * Plez (2012) Plez, B. 2012, Turbospectrum: Code for spectral synthesis
  * Ryabchikova et al. (2015) Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005
  * Sneden et al. (2008) Sneden, C., Cowan, J. J., & Gallino, R. 2008, ARA&A, 46, 241
  * Starkenburg et al. (2017) Starkenburg, E., Martin, N., Youakim, K., et al. 2017, Monthly Notices of the Royal Astronomical Society, 471, 2587–2604
  * Yoon et al. (2016) Yoon, J., Beers, T. C., Placco, V. M., et al. 2016, ApJ, 833, 20
  * Youakim et al. (2020) Youakim, K., Starkenburg, E., Martin, N. F., et al. 2020, MNRAS, 492, 4986

##  Appendix A Additional figures and data.

The figures in the appendix are described in Sections 2.3 and 3.2. Figures 12
and 13 describe the way we developed our method, whereas Figures 14 and 15
detail the dependence of our method on effective temperature and relative
carbon abundance. In Table 3 we provide the entire dataset that we used to
produce Figure 4, so that the interested readers can apply our method
themselves.

Figure 12: Contour plots where the third dimension are iso-flux lines, in
particular for the $fr\mathrm{{}_{CaHK/H\beta}}$ flux ratio. Those
iso-$fr\mathrm{{}_{CaHK/H\beta}}$ lines are temperature-sensitive up to about
$\mathrm{[Fe/H}]=-3.5$ dex, and then they become metallicity-sensitive, hence
we can use them to find a relation from which we can infer $\mathrm{[Fe/H]}$.
Because the metallicity sensitivity starts at $\mathrm{[Fe/H}]=-3.5$, we
cannot distinguish iron abundances below that threshold, but rather identify
whether they are below or above it. Table 3: Noise-free
$fr\mathrm{{}_{CaHK/H\beta}}$ and $fr\mathrm{{}_{G/CaNIR}}$ flux ratios and
their respective astrophysical parameters.

a/a | $fr\mathrm{{}_{G/CaNIR}}$ | $fr\mathrm{{}_{CaHK/H\beta}}$ | $\mathrm{T_{eff}}$ | $\mathrm{\log g}$ | $\mathrm{[Fe/H]}$ | $\mathrm{[C/Fe]}$
---|---|---|---|---|---|---
|  |  | (K) | (dex) | (dex) | (dex)
0 | 0.353852 | 0.056188 | 4800.0 | 5.0 | -0.5 | -0.5
1 | 0.361829 | 0.056728 | 4800.0 | 5.0 | -0.5 | -1.0
2 | 0.365205 | 0.056951 | 4800.0 | 5.0 | -0.5 | -1.5
3 | 0.337023 | 0.054872 | 4800.0 | 5.0 | -0.5 | 0.0
4 | 0.288912 | 0.050126 | 4800.0 | 5.0 | -0.5 | 0.5
5 | 0.187117 | 0.038496 | 4800.0 | 5.0 | -0.5 | 1.0
6 | 0.122764 | 0.032674 | 4800.0 | 5.0 | -0.5 | 1.5
7 | 0.086716 | 0.033368 | 4800.0 | 5.0 | -0.5 | 2.0
8 | 0.324806 | 0.050360 | 4800.0 | 5.0 | 0.0 | -0.5
9 | 0.333351 | 0.050936 | 4800.0 | 5.0 | 0.0 | -1.0
10 | 0.336977 | 0.051179 | 4800.0 | 5.0 | 0.0 | -1.5
11 | 0.305943 | 0.048967 | 4800.0 | 5.0 | 0.0 | 0.0
12 | 0.227209 | 0.042632 | 4800.0 | 5.0 | 0.0 | 0.5
13 | 0.144015 | 0.034670 | 4800.0 | 5.0 | 0.0 | 1.0
14 | 0.097575 | 0.032297 | 4800.0 | 5.0 | 0.0 | 1.5
15 | 0.073156 | 0.034265 | 4800.0 | 5.0 | 0.0 | 2.0
… | … | … | … | … | … | …
8963 | 0.471143 | 0.097089 | 5000.0 | 1.5 | -7.0 | 0.0
8964 | 0.471170 | 0.097088 | 5000.0 | 1.5 | -7.0 | 0.5
8965 | 0.471178 | 0.097102 | 5000.0 | 1.5 | -7.0 | 1.0
8966 | 0.471137 | 0.097087 | 5000.0 | 1.5 | -7.0 | 1.5
8967 | 0.471192 | 0.097109 | 5000.0 | 1.5 | -7.0 | 2.0
8968 | 0.471005 | 0.097071 | 5000.0 | 1.5 | -7.0 | 2.5
8969 | 0.470596 | 0.097043 | 5000.0 | 1.5 | -7.0 | 3.0
8970 | 0.469710 | 0.096964 | 5000.0 | 1.5 | -7.0 | 3.5
8971 | 0.466643 | 0.096713 | 5000.0 | 1.5 | -7.0 | 4.0
8972 | 0.458394 | 0.095966 | 5000.0 | 1.5 | -7.0 | 4.5
8973 | 0.439977 | 0.094247 | 5000.0 | 1.5 | -7.0 | 5.0
8974 | 0.409192 | 0.090923 | 5000.0 | 1.5 | -7.0 | 5.5
8975 | 0.368257 | 0.084704 | 5000.0 | 1.5 | -7.0 | 6.0
333Figure 4, as well as our model, was produced from this dataset.

Figure 13: We fit Equation 2 in the metallicity-$fr\mathrm{{}_{CaHK/H\beta}}$
space for each temperature-surface gravity pair from our parameter space (see
Figure 2). Some of those fits are shown in the above panels.

Figure 14: Temperature dependence of the performance of our method is twofold.
When $\mathrm{[C/Fe]}\leq 0$, the uncertainty is slightly better for lower
temperatures up to $\mathrm{[Fe/H]}\sim-3.5$, and above that, it becomes
almost independent of $\mathrm{T_{eff}}$. When $\mathrm{[C/Fe]}>0$, the
pattern is the same as above for $\mathrm{[Fe/H]}$ up to $\sim-3.5$; whereas,
above -3 dex, lower temperatures render higher uncertainties.

Figure 15: Computation of $\sigma_{\mathrm{{[Fe/H]}}}$ for all temperatures
corresponding to $\mathrm{\log g}=3.5$ dex, and all bins of $\mathrm{[C/Fe]}$.
The manner in which the performance of our method depends on the relative
carbon abundances relates to whether or not $\mathrm{[C/Fe]\lessgtr 0}$ dex:
at and below a Solar value, $\sigma_{\mathrm{{[Fe/H]}}}$ is practically
independent of carbon, but above it the uncertainty rises as carbon-
enhancement increases. The last effect lessens as temperature rises. The same
exercise was performed for all surface gravities of our parameter space, and
the results were similar.


# Pulse shape discrimination using a convolutional neural network for organic
liquid scintillator signals

K. Y. Jung  B. Y. Han  E. J. Jeon  Y. Jeong  H. S. Jo  J. Y. Kim  J. G. Kim
Y. D. Kim  Y. J. Ko  M. H. Lee  J. Lee  C. S. Moon  Y. M. Oh11footnotetext:
corresponding author H. K. Park  S. H. Seo  D. W. Seol  K. Siyeon  G. M. Sun
Y. S. Yoon  I. Yu

###### Abstract

A convolutional neural network (CNN) architecture is developed to improve the
pulse shape discrimination (PSD) power of the gadolinium-loaded organic liquid
scintillation detector to reduce the fast neutron background in the inverse
beta decay candidate events of the NEOS-II data. A power spectrum of an event
is constructed using a fast Fourier transform of the time domain raw waveforms
and put into CNN. An early data set is evaluated by CNN after it is trained
using low energy $\beta$ and $\alpha$ events. The signal-to-background ratio
averaged over 1-10 MeV visible energy range is enhanced by more than 20% in
the result of the CNN method compared to that of an existing conventional PSD
method, and the improvement is even higher in the low energy region.

##  1 Introduction

The signature of a prompt positron and a delayed neutron capture signal pair
from the inverse beta decay (IBD) reaction [1] enables us to detect the
electron antineutrino’s weak interaction signal buried in a large number of
background events. One of the major backgrounds which mimic the IBD signal
pair in the organic scintillation detector is the fast neutron’s moderation
and capture. In the case of the reactor experiments, it is preferred to deploy
the detector underground to reduce the cosmogenic neutron rate. However, very
short baseline experiments to solve the so-called Reactor Antineutrino Anomaly
(RAA) [2] must deploy their detectors at distances within tens of meters from
the reactor cores, where only a minimal shield from the cosmogenic background
is provided. Although most cosmogenic fast neutron background, which is not
related to the reactor operation, can be eliminated by subtracting the
reactor-off background data, there may remain uncertainty from the variation
of the fast neutron flux depending on the atmospheric conditions [3, 4]. To
reduce the systematic uncertainty that embeds this time-varying background
improve statistical precision, it is better to reduce the fast neutron
background as much as possible. With the help of digital signal processing,
the fast-$n$ background can be mitigated using the pulse shape discrimination
(PSD) at the analysis level. PSD is based on the fact that the scintillation
pulse shape of the fast-$n$’s moderation, mainly by proton recoil, is
different from that of $\gamma$’s or $e^{+/-}$’s electron recoil signal. Since
the difference in the pulse shapes lies primarily in the slow components of
the scintillation, the ratio of the pulse tail area to the total area
($F_{\textrm{tail}}$) in the time domain is generally used as the PSD
parameter. The discrimination power worsens at the lower visible energy side
due to poorer photostatistics and signal-to-noise ratio. Other methods using
the Fourier transform have been tested in order to get over the degradation by
high-frequency noises [5, 6, 7].

While the traditional PSD techniques mainly focus on finding principle
components and constructing simple parameters, convolutional neural networks
(CNN) can take the raw signals as the input data. For decades, CNN
architectures have been introduced and developed for image classification [8,
9, 10, 11, 12]. In the neutrino physics field, CNN has been studied for the
classification of neutrino interaction types, particle identification,
background rejection, and so on [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24].

NEOS [25] used a linear-alkylbenzene (LAB) based gadolinium (Gd)-loaded liquid
scintillator, into which a 10% volume of a di-isopropylnaphthalene (DIN) based
liquid scintillator was mixed to enhance the PSD capability [26, 27]. As a
result of the phase-I measurement using the $F_{\textrm{tail}}$ method, the
background was reduced by 70%, and a signal-to-background ratio of 22 was
achieved in the 1-10 MeV prompt energy range [28, 25]. We have reported a
possible PSD improvement using a CNN applied to a test data set consisting of
fast-$n$ and $\beta/\gamma$ background events [29]. The time domain waveforms
used to calculate $F_{\textrm{tail}}$ were given as input data for the CNN
training and evaluation.

We have continued the study for further PSD improvement and to see the impact
on the experimental result of NEOS-II, which is presented in the following
sections. First, data processing for the conventional PSD method utilizing the
tail charge and the construction of the frequency domain power spectrum using
fast Fourier transform (FFT) in the NEOS-II analysis are described in detail
in section 2. In section 3 the CNN method using the FFT power spectrum is
explained. We will apply both the conventional and CNN methods on the NEOS-II
IBD candidate dataset and compare the results in section 4.

##  2 Data preparation for PSD

###  2.1 A conventional method

A signal in the liquid scintillator is detected by 2$\times$19 photomultiplier
tubes (PMTs) in the NEOS detector [25]. Each PMT signal is digitized via a
flash analog-to-digital converter (FADC) module which has a 12-bit resolution
for 2.5 volts peak-to-peak dynamic range and a 500 MHz sampling rate. The
digitized 38 waveforms in a 480 nanosecond window are then stored as the raw
data of an event. An example of the raw waveforms of an event is shown in
figure 1a.

For the pulse shape analysis in the time domain, a synchronized waveform
($w_{s}$) is constructed (see figure 1b). Considering the detector size and
the PMT transit time spread, signal timings in channels can be different by
several nanoseconds. First, we look at the rising part of the $i$th PMT’s raw
waveform between the maximum point and the last point exceeding the pedestal
root mean square before the maximum. When the rising part consists of more
than 3 points, we define a pulse time ($t_{p}^{i}$) of the waveform at the
half maximum point found by fitting
$\left(t,\,\left[-\log\left(w_{i}(t)/w_{i}^{\mathrm{max}}\right)\right]^{1/2}\right)$
to a linear function, where $w_{i}(t)$ stands for the waveform amplitude at a
time bin $t$ for $i$th PMT, and $w_{i}^{\mathrm{max}}$ for its maximum
amplitude. While the total charge of an event ($q_{\mathrm{sum}}$), a measure
of the energy deposit in the detector [28], is calculated from the sum of the
integrated area of all 38 raw waveforms, a synchronized waveform ($w_{s}$) of
the event is formed by accumulating the waveforms of which $t_{p}$’s are
defined:

$w_{s}(t)=\sum_{i}w_{i}\left(t_{p}^{i}-25+t\right),$ (2.1)

where $t$ is the time bin index from 1 to 170. The selection of 170-time bins
allows it to contain most of the main signal information and exclude the
pedestal and possible afterpulses. The tail-to-total charge ratio,
$F_{\mathrm{tail}}$, is calculated from $w_{s}$. In the NEOS-II analysis, we
have developed another PSD parameter, which is defined as the spread of mean-
time in the tail range of $w_{s}(t)$ as follows:

$\sigma_{\mathrm{tail}}(t_{s})=\sqrt{\sum_{t=t_{s}}^{170}\frac{t^{2}\cdot
w_{s}(t)}{A_{\mathrm{total}}}-\sum_{t=t_{s}}^{170}\left[\frac{t\cdot
w_{s}(t)}{A_{\mathrm{total}}}\right]^{2}},$ (2.2)

where $t_{s}$ is the tail starting bin, and
$A_{\mathrm{total}}=\sum_{t=1}^{170}w_{s}(t)$ is the total area of the
synchronized waveform.

Figure 1: An example of the raw data processing of an event for PSD. (a) Raw
waveforms ($w_{i}$’s) from 38 PMTs. All $x$-axes share the same time range
(0-480 ns), and a division in $y$-axes stands for 5 mV. The pulse time
($t_{p}$), drawn as an orange cross, is defined at the half-maximum point in
the rising part. (b) Waveforms whose $t_{p}$’s are properly found are
synchronized and accumulated to produce a synchronized waveform (blue).
$F_{\mathrm{tail}}$ (orange) and $\sigma_{\mathrm{tail}}$ (green) are then
calculated for each tail-starting time bin. (c) Power spectra of 38 PMT
waveforms (thin gray) and their sum (thick black). See the text for a detailed
description.

The conventional optimization of PSD is finding $t_{s}$ which maximizes the
discrimination power, i.e., the figure of merit (FoM), between $\beta/\gamma$
and fast-neutron/$\alpha$ events:

$\mathrm{FoM}=\frac{\left|m_{\gamma}-m_{n}\right|}{\sqrt{\sigma_{\gamma}^{2}+\sigma_{n}^{2}}},$
(2.3)

where $m$ and $\sigma$ denote the Gaussian mean and width of a type of
scintillation event–$\gamma$ for $\beta/\gamma$ and $n$ for fast-
neutron/$\alpha$, respectively. Figure 2a shows how FoM and optimum $t_{s}$
values change for different $q_{\mathrm{sum}}$ ranges.

Figure 2: (a) Discrimination powers (FoM) at different tail starting points
for different charge sum ranges; (b) PSD parameter, $\sigma(t_{s}=46)$ versus
charge sum; (c) $\sigma(t_{s}=46)$ distribution and its fit to a two-Gaussian
function for events in 2-3 MeV range (charge sum of 1000-1500 pC range). Data
shown here are the prompt events, followed by $n$-capture-on-Gd events in the
reactor-off period data.

It is natural to have a small FoM in the low $q_{\mathrm{sum}}$ range, as
shown in figure 2, because of a small number of photons, which may be further
degraded by choice of waveforms with valid $t_{p}$’s. An example of a low
energy ($\sim$1 MeV, or charge 500 pC) event in figure 1 shows that $t_{p}$
cannot be found for some channels (10, 11, 34, and 35).

###  2.2 Utilization of the fast Fourier transform

Generally, a proton recoil event shows a larger decay time in the slow
decaying component in the organic scintillator than an electron recoil event.
It means that, when the total sizes of the two pulses are the same, the pulse
of the former varies smaller in the decaying part and will show a smaller
amplitude at the corresponding frequency range after the Fourier transform. A
simple example is shown in figure 3.

Figure 3: Two exponential functions with different decay times (left) and
their normalized squared amplitudes in the frequency domain after a fast
Fourier transform (right). The low frequency region is shown on a linear scale
in the inset.

A normalized power spectrum in the frequency domain, denoted as
$\mathcal{P}(\omega)$, can be constructed by adding the amplitude-squared of
the Fourier transform ($\mathcal{F}$) of each PMT’s waveform:

$\mathcal{P}(\omega)=\sum_{i=1}^{38}\left|\mathcal{F}\left(w_{i}(t)\right)\right|^{2}.$
(2.4)

The fast Fourier transform is used since the waveforms have a finite length
and are discretized in time. The waveforms in the time bins between the 25th
from the start and the 20th from the end are used for the transform, to
exclude useless pedestal parts and to avoid possible afterpulses,
respectively. Figure 1c shows each of the 38 power spectra
($|\mathcal{F}(w_{i}(t))|^{2}$) and their sum ($\mathcal{P}(\omega)$). The
zeroth term of $\mathcal{P}(\omega)$, $\mathcal{P}(0)$, is proportional to the
quadratic sum of the integrated area of the time domain signal and can be used
as the denominator for normalization. One advantage of using the frequency
domain power spectrum is that one doesn’t have to worry about time
information, i.e., finding pulse times that may introduce an additional error
so that all 38 PMT signals can be used. Another advantage is data reduction.
Neglecting the phase terms leaves a half of the number of original time bins.
It helps save computing power when the spectrum is used as an input for CNN.

###  2.3 Preparation of data samples

There are many ways to get pure electron recoil events in the NEOS data, while
it is almost impossible to sort out pure proton recoil events in the wide
energy range of interest ($\sim$1$<E$ /MeV$<$10). Gamma-event dominant data
can be found in most of the source calibration data. Neutron-capture $\gamma$
events depositing energies up to around 8 MeV in the detector can be found in
the delayed coincident event(s) of 252Cf or PoBe source calibration data. In
the case of 252Cf, multiple neutrons’ moderation and capture may happen in an
event time window, so that pile-up signals can pollute some waveforms. In the
PoBe case, the source emits a fast neutron with or without a 4.4 MeV $\gamma$.
A proton recoil by the fast neutron’s moderation leaves a visible energy of
less than 2 MeV, and sometimes a pile-up happens by the coincidence of the 4.4
MeV $\gamma$ and the moderation of the neutron. Another drawback is that
events mostly occur at around the top and middle of the horizontal cylindrical
detector due to the source position, so the scintillation photons are observed
not uniformly or symmetrically for 38 PMTs but mainly by PMTs at the top part.
Other calibration sources, 22Na, 60Co, and 137Cs, provide symmetric charge
distributions of $\gamma$ events for 38 PMTs but are limited to the low energy
side.

In the background data, events in the muon veto time window, 150 microseconds
after a muon counter event, have rich information, although they are
unselected in the IBD reconstruction. First, the Michel electron [30] can be
found by a delayed coincidence with stopping muon events. Some of them may
coincide with afterpulses of the large muon signals; the contribution of the
afterpulse becomes negligible at high energy. Therefore the Michel electron
event above 4 MeV can be used as an electron recoil sample. Second, there are
fast neutrons following a muon. We can find the moderation of fast neutrons
leaving proton recoil signals and the electron recoil signals of $\gamma$
events by the neutron captures (figure 2). In this case, both types of recoil
events leave visible energy wide enough for our interest, but the events
cannot be clearly classified into the recoil types without PSD.

Consecutive $\beta-\alpha$ decays of 214Bi and 214Po, originating from the
radon contamination in the liquid scintillator, provide clean samples of
electron and proton recoil events, respectively, at the low visible energy
range. Thanks to the short half-life of 214Po (164 microseconds) and the full
energy deposit from its $\alpha$ decay, the pair of decay events can be
selected by a 300 microsecond coincidence time window, a cut on the PMT charge
sum around 500 pC for the delayed event, and a loose $F_{\mathrm{tail}}$
selection, as shown in figure 5a. Then a fine selection of the charge and
energy values is made for the $\alpha$ events.

Table 1: Usable data samples for the PSD analysis. Calibration sources have
typical radioactivities of $\mathcal{O}(10^{2})$ Bq.

Data type | Source | Recoil | Visible energy (MeV) | Comment
---|---|---|---|---
Calibration |  22Na,60Co,137Cs | $e$ |  $<\sim$2.5 | Low energy
252Cf | $e$ |  $<$10 | Pile-up of multiple neutrons
$p$ | DR∗
PoBe | $e$ | 1-10 | Asymmetric source position
$p$ | <2 | Pile-up of prompt fast-$n+\gamma$
Background | Michel electron | $e$ |  $>$4 | High energy, $\sim$0.7 Hz
$\mu$-induced $n$ | $e$ | 1-10 | Similar to IBD
$p$ | DR∗ |  $\mathcal{O}(10^{2})$ events/day
Bi($\beta$)-Po($\alpha$) | $e$ |  $<$3 | Low energy, clean samples
$p$ |  $\sim$1 |  $\mathcal{O}(10^{4})$ events/day

∗DR: in the whole dynamic range of the detector: 0.6-20 MeV.

##  3 PSD using Convolutional Neural Network

###  3.1 General framework

A neural network optimizes weight parameters that transform input data within
the hidden layers of the network. The traditional deep learning connects the
weight parameters to the full input shape, which may cause overfitting and a
lack of generalization. The weight parameters in CNN are connected to a subset
of the input shape and shared by the whole input to avoid overfitting and to
improve generalization [31]. There is a series of hyperparameters when CNN
optimizes the weight parameters, for example, kernel, filter, stride,
activation function, pooling, dropout, etc. With these hyperparameters, CNN
extracts the $j$’th feature map, $\hat{x}^{l}_{j}$, defined as

$\hat{x}^{l}_{j}\equiv
f\left(\overset{m}{\underset{i=1}{\sum}}x_{i+j}w^{l}_{i}+b^{l}_{j}\right),$
(3.1)

where $m$ is a given kernel size, $l$ is the number of filters, $x_{i}$ is the
$i$’th input, $f$ is an activation function, b is a bias, and $w$ represents
weights. The kernel is the size of shared parameters in a layer. The weights
in the kernel are combined with input data in an activation function to
extract the feature map. A filter is a number of kernels and a depth of
feature map to learn various features. The stride indicates the amount of
movement for the kernel over the input. The activation function is an operator
on the linear combination of input and weight parameters to produce features;
ReLU (rectified linear unit), hyperbolic tangent, and sigmoid functions are
used in this study. Pooling and dropout are used to prevent overfitting which
is specialized for CNN in specific samples. Pooling is a way to reduce the
size of a feature map, such as max pooling, average pooling, and stochastic
pooling. In this study, max pooling is adopted, choosing a max value of the
neighbor in the feature map. Dropout removes a part of the feature map to
prevent overfitting and improve the generality of CNN. Figure 4 shows the
general framework of CNN and each hyperparameter’s role. When the CNN has more
than one layer, the feature map after the convolutional layer becomes the
input of the next layer.

For the supervised binary classification, each input getting into the CNN
training needs to be labeled as the true information. In our case, the label
for an electron recoil ($\beta/\gamma$) event is assigned as 1, and 0 for a
proton recoil ($\alpha/$fast-$n$) event. Each event is evaluated with a score
of a CNN result between 0 and 1, indicating how close the event is to one
label or the other. A loss is calculated from the difference between the score
and the label. The binary entropy function is defined as:

$L(w)=-\frac{1}{N}\sum_{j=1}^{N}\left[y_{j}\cdot
log(\hat{y}_{j})+(1-y_{j})\cdot log(1-\hat{y}_{j})\right],$ (3.2)

where $\hat{y}_{j}$ and $y_{j}$ are the scores as a function of weights and
the label for the $j$-th input. This loss function can diverge infinitely when
the difference between the label and score gets larger. CNN updates the
optimized weights in the direction of gradient descent of the loss function by
backpropagation.

Figure 4: An example of CNN framework and hyperparameters’ role. This cartoon
shows three kernels, four filters, max pooling of two, and a dropout of 0.5.
The index is the same as in equation 3.1

###  3.2 CNN training and architecture

It is our best choice to use the Bi-Po events described in section 2.3 as the
training input for CNN because they provide us with the cleanest samples of
the electron and proton recoil events at low energy, where an efficient
rejection of the fast-$n$ background in the IBD candidate events is difficult
by the conventional PSD method. Similarities of the synchronized waveforms
between $\beta$ and $\gamma$ and between $\alpha$ and fast-$n$ are already
shown in our previous work [29]. The normalized power spectrum,
$\mathcal{P}(\omega)/\mathcal{P}(0)$, which has 98 data points with its label,
is put into the CNN input layer. The reasons why the FFT spectrum is used
instead of the time-synchronized waveform are; (1) to fully utilize the 38 PMT
waveforms as explained in section 2.2, (2) to avoid the possible bias caused
by setting the pulse time at a given time point for the synchronized waveform,
and (3) to save computing resources by reducing the input data shape from 170
of the time synchronized waveform to 98 of the FFT power spectrum. Figure 5a
shows the selected $\beta$ and $\alpha$ events in the $F_{\mathrm{tail}}$
versus charge plane, and figure 5b shows the averages of their normalized
power spectra. One can see that the main difference between the two power
spectra lies in the low-frequency range.

Considering separation power and generality, we tried constructing CNN
architectures as simply as possible using large strides and a few filters
[32]. Because there is no mathematical formula to find an optimized CNN
architecture, the appropriate architecture and hyperparameters are found by
trial and error [33]. Out of many different CNN architectures that are tried,
the ones with relatively simple layers and structures for 1D CNN show
satisfactory performances in signal classification compared to more
complicated architectures in terms of accuracy and loss. One of the simplest
architectures that offers the lowest loss and the best accuracy for the
evaluation data samples is chosen.

Figure 5: Consecutive 214Bi’s $\beta$\- and 214Po’s $\alpha$\- decay event
candidates. (a) $F_{\mathrm{tail}}$ versus charge distribution. The $\beta$
events are selected for $F_{\mathrm{tail}}<0.16$, and the $\alpha$ events are
selected inside the black dashed ellipse. (b) Averages (solid curves) and the
standard deviations of normalized power spectra for $\beta$ (blue) and
$\alpha$ (orange) events where $\beta$ events are selected to have a similar
mean charge with the $\alpha$ events. Figure 6: The CNN architecture design
for this work. See the text for the detailed description.

The CNN architecture of this work is drawn in figure 6. Keras [34] and
Tensorflow [35] are used for CNN training on a computer with an NVIDIA GeForce
RTX 2080 GPU. Data for training are divided into 80% for the training set and
20% for the validation set. First, every batch consisting of 128 events, each
of which has 98 data points, in the training data set goes through the first
convolutional layer consisting of 16 kernels of size three and stride 3. Its
outputs are then put into the second convolutional layer with 32 kernels of
size two and stride 2. During these feature extraction processes, the ReLU
activation function is applied to every convolutional layer. A max-pooling and
a dropout are applied to the second layer and every convolutional layer,
respectively, to prevent overfitting. The feature maps extracted through the
convolutional and pooling layers are then fully connected to 128 weights by a
$\tanh$ activation function for the feature classification. Finally, the
result from the fully-connected layer turns into the score by the sigmoid
function. The weights are updated in the direction of loss minimization for
the next batch, using the binary cross entropy (equation 3.2) and an Adam
optimizer [36]. An epoch is completed when the entire batches pass through the
CNN training, where computing time for each epoch takes less than 30 sec.

###  3.3 Training result

Overfitting may occur when the loss of the training set decreases while the
loss of the validation grows. It is also considered overfitting and the
training is stopped when the validation loss does not decrease within five
epochs. The optimized weights are chosen at the epoch right before the
overfitting happens. We expect the overfitting to occur by the fortieth epoch
and set the hard limit there. As shown in figure 7, the training does not
improve the validation loss after the seventh epoch. As a result, we obtained
99.93% accuracy for the validation data set at epoch 7.

Figure 7: (a) Accuracy (black) or loss (red) during the epoch. Solid and dot
lines represent the training set and validation set, respectively. CNN model
is selected at epoch 7 (orange dashed line) to prevent overfitting. (b) Score
distribution in the validation set. Bins are spaced using the error function:
$0.5\cdot[1+\mathrm{erf}\left((x-0.5)/0.18\right)]$ for equally spaced $x$
between 0 and 1, to effectively show the dense distribution at each end.

##  4 Results

###  4.1 Score cut and efficiency

With the conventional PSD method, our selection strategy for the IBD candidate
events is to maintain the selection efficiency for the electron recoil events
at a value higher than 99% and to have moderate rejection efficiencies for the
fast-$n$ background throughout the energy range of interest (1-10 MeV). It is
to minimize any distortion and systematic uncertainty in the IBD prompt energy
spectrum caused by the selection efficiency. This strategy, however, may
introduce an additional systematic error for the unknown difference between
the measured background from reactor-off periods and the background mixed with
IBD in the reactor-on data. Remained fast-$n$ background rate by an
inefficient rejection, especially at low energy, can change in time due to an
unstable detector response [37] or by fluctuations of the environmental
condition [4]. In the NEOS-I case with the conventional PSD method, a 10%
conservative error for the reactor-off spectrum was assigned due to this
unknown background fluctuation [25].

Our goal in employing CNN for IBD selection is to maintain the same high
selection efficiency for the electron recoil events while improving the
fast-$n$ rejection efficiency. The selection efficiency is investigated using
evaluated CNN scores of the clean electron recoil events from the 22Na source
calibration data and the Michel electron data for the low and high energy
ranges, respectively. In this work, the data are sorted into different energy
windows and a score of 99.5% acceptance for each energy window is found, as
shown in figure 8. A score cut as a function of energy is found by
interpolation of those score points. The cut value gets close to 1 as the
energy increases.

Figure 8: CNN score distributions of 22Na source data (top) and Michel
electron events (bottom) at different energy ranges. Binning is the same as in
figure 7b. The black dashed lines and the numbers indicate the score cuts
satisfying 99.5% $\gamma$ selection efficiency.

As mentioned in section 2.3, clean fast-$n$-only events, especially at low
energy, cannot be sorted out in the NEOS data. Estimation of efficiency with
the CNN score is not as easy for the fast-$n$ rejection as for the
$\beta$/$\gamma$ selection. Therefore, the CNN score cut in this work is
validated by checking consistency among a result by CNN, another by
conventional PSD, and the other without PSD. Since the $\beta$/$\gamma$
selection efficiencies were set to be higher than 99% in both the CNN and the
conventional methods, the reactor on-off spectra from all three methods should
agree.

###  4.2 Application to the IBD selection

PSD parameter ($\sigma_{\mathrm{tail}}$) of the IBD prompt event candidates
from an early reactor-on data set is shown in figure 9 The
$\sigma_{\mathrm{tail}}$ of $\beta$/$\gamma$ events lies along around 16, and
the resolution becomes poorer as energy decreases. The
$\sigma_{\mathrm{tail}}$ cut value along with the fast-$n$’s
$\sigma_{\mathrm{tail}}$ decreases as energy grows. The cut for this early
data set seems quite loose. Still, it is to sustain the same efficiency with
the following data sets with poorer discrimination caused by a light yield
degradation during the NEOS-II measurement [37]. In terms of CNN score, most
events are separated into scores very close to 0 or 1, and not many events
exist in the middle. Some events on the high-energy side and above the
$\sigma_{\mathrm{tail}}$ cut get scores close to half, mainly because the CNN
is optimized using only the low-energy Bi-Po events. However, events with dark
red color can only survive after the CNN cut for the high-energy side because
the cut is at a score higher than 0.99 for $E>4$ MeV, as shown in figure 8. On
the low energy side, in contrast, a considerable number of events below the
$\sigma_{\mathrm{tail}}$ cut are classified clearly as fast-$n$ events. In
short, the CNN score cut rejects more events, many of which are in the low-
energy region.

Figure 9: The left figure shows the PSD parameter ($\sigma_{\mathrm{tail}}$)
versus the energy of the IBD prompt event candidates in the early reactor-on
data of NEOS-II. Color represents the evaluated CNN score. The magenta dashed
curve indicates a PSD cut on the conventional PSD parameter
($\sigma_{\mathrm{tail}}$). The right figure shows the projected data on
$\sigma_{\mathrm{tail}}$. Blue and orange histograms of solid (dashed) lines
are selected $\beta$/$\gamma$ and fast-$n$ events, respectively, by the CNN
score ($\sigma_{\mathrm{tail}}$) cut.

Figure 10 shows the IBD prompt energy distributions of the signals (reactor-on
minus reactor-off) and backgrounds (reactor-off) with three analysis methods:
one without PSD, another with $\sigma_{\mathrm{tail}}$ cut, and the other by
the CNN score cut. Both PSD methods show significant backgroud reductions,
while the CNN method rejects more background at low energy. Despite
significant differences among the background spectra, subtracting them from
their corresponding reactor-on spectra leaves the IBD signal spectra with
almost perfect agreements among one another. The total signal rate (on$-$off)
after the $\sigma_{\mathrm{tail}}$ cut is 100%, and the one by CNN is 99.8%,
compared to the result of 1864 events per day without PSD, respectively. Daily
background rates from CNN and $\sigma_{\mathrm{tail}}$ methods are
59.9$\pm$0.7 and 76.2$\pm$0.8, respectively, enhancing the signal-to-
background ratio by more than 20%, from 24 to 31. As intended and expected,
improvement is more significant on the low-energy side as CNN is trained with
low-energy data.

Figure 10: Top: IBD signal (ON-OFF; solid line) and background (OFF; dashed
line) spectra with different analysis methods–no PSD (gray), conventional PSD
(green), and PSD using CNN (blue). Bottom: Signal-to-background ratio for each
analysis method.

##  5 Summary and discussion

This study shows that CNN helps PSD to improve fast-$n$ background rejection
by a considerable amount, at least in the case that a suitable data set exists
for training in the region where the conventional PSD shows limited power.
Excellent separation between $\beta$/$\gamma$ and fast-$n$/$\alpha$ events may
let us avoid complicated cut design and efficiency estimation in the
conventional method. At the same time, it takes a similar time and effort to
optimize the CNN architecture. It is also shown that a normalized raw FFT
power spectrum can be a good input for CNN without much loss of key features
of recoil event type and avoiding bias.

Figure 11: The figure of merit (FoM) for the PSD parameter for the prompt
events of 2-4 MeV energy during the reactor-on period of the NEOS-II data
taking. The decrease in FoM is due to the continuous light yield degradation.
The data are divided into 9 groups and the first data group of 44 days and
preceding reactor-off data for 29 days are shown in this study.

The result presented here is for the first 29+44 days with the reactor off+on,
respectively, out of the total 500 days of NEOS-II data taking. The FoM of the
PSD parameter became lower in time as shown in figure 11, because of the
degradation of the scintillation light yield during the whole NEOS-II data-
taking period, as reported in [37]. Therefore, data are divided into 9 groups
for the reactor-on period, and each data group is trained and evaluated
separately using the same process developed in this study. A fine-tuning of
the CNN score cut for each data group is going on to keep the selection
efficiency to be as stable and close to unity as possible to minimize the
uncertainty from the efficiency fluctuation. A preliminary result for the
whole NEOS-II data set has been presented in [37]. One of the primary goals of
NEOS-II is to differentiate the IBD spectrum of 235U and that of plutonium.
This spectrum decomposition is derived from the change of the total spectrum
in time, which is only a few percent for a reactor operation cycle. Therefore,
small improvements in statistical and systematic accuracy will help improve
the experimental result more than a little. However, the signal-to-background
ratio is already as large as 24 without the help of CNN.

A limited feature of this study is that only low visible energy $\alpha$
events are used in the CNN training to mitigate the fast-$n$ background whose
visible energy range spans all the IBD prompt energy ranges and above. Further
PSD improvement is expected to be achieved, at least by getting clean fast-$n$
and $\beta$/$\gamma$ data samples at the rest of the energy range.

## Acknowledgments

This work is supported by the National Research Foundation of Korea
(NRF-2022R1A2C1009686). NEOS-II is supported by Institute for Basic Science
(IBS-R016-D1 and IBS-R016-D1-2018-B01). Computational works for the research
were performed on the data analysis hub, Olaf in the IBS Research Solution
Center.

## References

  * [1] P. Vogel and J.F. Beacom, _Angular distribution of neutron inverse beta decay, anti-neutrino(e) + p — $>$ e+ + n_, _Phys. Rev. D_ 60 (1999) 053003 [hep-ph/9903554].
  * [2] G. Mention, M. Fechner, T. Lasserre, T.A. Mueller, D. Lhuillier, M. Cribier et al., _Reactor antineutrino anomaly_ , _Physical Review D_ 83 (2011) .
  * [3] J.R. Davis, E. Brubaker and K. Vetter, _Fast neutron background characterization with the Radiological Multi-sensor Analysis Platform (RadMAP)_ , _Nucl. Instrum. Meth. A_ 858 (2017) 106 [1611.04996].
  * [4] STEREO collaboration, _Improved sterile neutrino constraints from the STEREO experiment with 179 days of reactor-on data_ , _Phys. Rev. D_ 102 (2020) 052002 [1912.06582].
  * [5] G. Liu, M.J. Joyce, X. Ma and M.D. Aspinall, _A digital method for the discrimination of neutrons and $\gamma$ rays with organic scintillation detectors using frequency gradient analysis_, _IEEE Transactions on Nuclear Science_ 57 (2010) 1682.
  * [6] M.J. Balmer, K.A. Gamage and G.C. Taylor, _Comparative analysis of pulse shape discrimination methods in a 6li loaded plastic scintillator_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 788 (2015) 146.
  * [7] M. Hubbard, M. Taggart and P. Sellin, _Exploration of Fourier based algorithms and detector designs for pulse shape discrimination_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 930 (2019) 64.
  * [8] Y. Lecun, L. Bottou, Y. Bengio and P. Haffner, _Gradient-based learning applied to document recognition_ , _Proceedings of the IEEE_ 86 (1998) 2278.
  * [9] A. Krizhevsky, I. Sutskever and G.E. Hinton, _Imagenet classification with deep convolutional neural networks_ , _Advances in neural information processing systems_ 25 (2012).
  * [10] K. Simonyan and A. Zisserman, _Very deep convolutional networks for large-scale image recognition_ , 2014.  10.48550/arxiv.1409.1556.
  * [11] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov et al., _Going deeper with convolutions_ , 2014.  10.48550/ARXIV.1409.4842.
  * [12] K. He, X. Zhang, S. Ren and J. Sun, _Deep residual learning for image recognition_ , 2015.  10.48550/ARXIV.1512.03385.
  * [13] A. Aurisano, A. Radovic, D. Rocco, A. Himmel, M.D. Messier, E. Niner et al., _A Convolutional Neural Network Neutrino Event Classifier_ , _JINST_ 11 (2016) P09001 [1604.01444].
  * [14] DUNE collaboration, _Neutrino interaction classification with a convolutional neural network in the DUNE far detector_ , _Phys. Rev. D_ 102 (2020) 092003 [2006.15052].
  * [15] MicroBooNE collaboration, _Convolutional Neural Networks Applied to Neutrino Events in a Liquid Argon Time Projection Chamber_ , _JINST_ 12 (2017) P03011 [1611.05531].
  * [16] MicroBooNE collaboration, _Deep neural network for pixel-level electromagnetic particle identification in the MicroBooNE liquid argon time projection chamber_ , _Phys. Rev. D_ 99 (2019) 092001 [1808.07269].
  * [17] MicroBooNE collaboration, _Semantic segmentation with a sparse convolutional neural network for event reconstruction in MicroBooNE_ , _Phys. Rev. D_ 103 (2021) 052012 [2012.08513].
  * [18] NEXT collaboration, _Background rejection in NEXT using deep neural networks_ , _JINST_ 12 (2017) T01004 [1609.06202].
  * [19] N. Choma, F. Monti, L. Gerhardt, T. Palczewski, Z. Ronaghi, Prabhat et al., _Graph neural networks for IceCube signal classification_ , _CoRR_ abs/1809.06166 (2018) [1809.06166].
  * [20] E. Racah, S. Ko, P. Sadowski, W. Bhimji, C. Tull, S.-Y. Oh et al., _Revealing fundamental physics from the Daya Bay neutrino experiment using deep neural networks_ , in _2016 15th IEEE International Conference on Machine Learning and Applications (ICMLA)_ , pp. 892–897, 2016, DOI.
  * [21] DeepLearnPhysics collaboration, _Scalable deep convolutional neural networks for sparse, locally dense liquid argon time projection chamber data_ , _Phys. Rev. D_ 102 (2020) 012005 [1903.05663].
  * [22] KamLAND-Zen collaboration, _Search for Majorana Neutrinos near the Inverted Mass Hierarchy Region with KamLAND-Zen_ , _Phys. Rev. Lett._ 117 (2016) 082503 [1605.02889].
  * [23] J. Griffiths, S. Kleinegesse, D. Saunders, R. Taylor and A. Vacheret, _Pulse Shape Discrimination and Exploration of Scintillation Signals Using Convolutional Neural Networks_ , 1807.06853.
  * [24] MINERvA collaboration, _Reducing model bias in a deep learning classifier using domain adversarial neural networks in the MINERvA experiment_ , _JINST_ 13 (2018) P11020 [1808.08332].
  * [25] NEOS collaboration, _Sterile Neutrino Search at the NEOS Experiment_ , _Phys. Rev. Lett._ 118 (2017) 121802 [1610.05134].
  * [26] B.R. Kim, B.Y. Han, E.J. Jeon, K.K. Joo, J. Kang, N. Khan et al., _Pulse shape discrimination capability of metal-loaded organic liquid scintillators for a short-baseline reactor neutrino experiment_ , _Physica Scripta_ 90 (2015) 055302.
  * [27] NEOS collaboration, _Development and Mass Production of a Mixture of LAB- and DIN-based Gadolinium-loaded Liquid Scintillator for the NEOS Short-baseline Neutrino Experiment_ , _J. Radioanal. Nucl. Chem._ 310 (2016) 311 [1511.05551].
  * [28] NEOS collaboration, Y.M. Oh et al., _NEOS result and prospects_ , June 2018.  10.5281/zenodo.1286994.
  * [29] NEOS II collaboration, _Pulse-shape Discrimination of Fast Neutron Background using Convolutional Neural Network for NEOS II_ , _J. Korean Phys. Soc._ 77 (2020) 1118 [2009.13355].
  * [30] L. Michel, _Interaction between four half spin particles and the decay of the $\mu$ meson_, _Proc. Phys. Soc. A_ 63 (1950) 514.
  * [31] I. Goodfellow, Y. Bengio and A. Courville, _Deep Learning_ , MIT Press (2016).
  * [32] A. Khan, A. Sohail, U. Zahoora and A. Saeed, _A survey of the recent architectures of deep convolutional neural networks_ , _Artificial Intelligence Review_ 53 (2020).
  * [33] J. Rala Cordeiro, A. Raimundo, O. Postolache and P. Sebastião, _Neural architecture search for 1d cnns—different approaches tests and measurements_ , _Sensors_ 21 (2021) 7990.
  * [34] F. Chollet et al., _Keras_ , 2015.
  * [35] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean et al., _Tensorflow: A system for large-scale machine learning_ , _CoRR_ abs/1605.08695 (2016) [1605.08695].
  * [36] D.P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , 2014\.  10.48550/ARXIV.1412.6980.
  * [37] J. Kim, _Sterile neutrino i_neos-ii new results_ , June 2022.  10.5281/zenodo.6680618.