XXXX-XXXX

1]Department of Physics, Okayama University, Okayama, Japan 2]Department of
Physics and Astronomy, University of Utah, Salt Lake City, USA 3]Department of
Physics and Astronomy, University of Padova, Padova, Italy 4]Ecole
Polytechnique, Paris, France 5]National Astronomical Observatory of Japan,
Mitaka, Japan 6]TRIUMF, Vancouver, Canada 7]Research Center for Nuclear
Physics, Osaka, Japan 8]Department of Physics, Kyoto University, Kyoto, Japan

# Measurement of $\gamma$-rays generated by neutron interaction with 16O at 30
MeV and 250 MeV

T. Tano  T. Horai  Y. Ashida  Y. Hino  F. Iacob  A. Maurel  M. Mori  G.
Collazuol  A. Konaka  Y. Koshio  T. Nakaya  T. Shima  R. Wendell  [  [  [  [
[  [  [  [<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Deep understanding of $\gamma$-ray production from the fast neutron reaction
in water is crucial for various physics studies at large-scale water Cherenkov
detectors. We performed test experiments using quasi-mono energetic neutron
beams ($E_{n}=30$ and 250 MeV) at Osaka University’s Research Center for
Nuclear Physics to measure $\gamma$-rays originating from the neutron-oxygen
reaction with a high-purity germanium detector. Multiple $\gamma$-ray peaks
which are expected to be from excited nuclei after the neutron-oxygen reaction
were successfully observed. We measured the neutron beam flux using an organic
liquid scintillator for the cross section measurement. With a spectral fitting
analysis based on the tailored $\gamma$-ray signal and background templates,
we measured cross sections for each observed $\gamma$-ray component. The
results will be useful to validate neutron models employed in the on-going and
future water Cherenkov experiments.

xxxx, xxx

##  1 Introduction

A precise modeling of the neutrino-oxygen neutral-current quasi-elastic (NCQE)
interactions is essential for various studies of astrophysics and particle
physics at large-scale water Cherenkov detectors, such as Super-Kamiokande
(SK) bib:superk  and its successor, Hyper-Kamiokande bib:hyperk . For example,
atmospheric neutrino NCQE interactions on ${\rm{}^{16}O}$ in water needs to be
estimated accurately as they mimic the inverse beta decay signal
($\bar{\nu}_{e}+p\to e^{+}+n$) of the diffuse supernova neutrino background
(DSNB) 2021PhRvD.104l2002A ; Super-Kamiokande:2023xup . There have been
multiple neutrino-oxygen NCQE cross section measurements performed by using
accelerator neutrino beams in T2K 2014PhRvD..90g2012A ; 2019PhRvD.100k2009A
and atmospheric neutrinos at SK 2019PhRvD..99c2005W ; 2024PhRvD.109a1101S . At
neutrino energies of $E_{\nu}\gtrsim 200$ MeV, a single nucleon knock-out is
expected to occur dominantly after the NCQE interaction 2012PhRvL.108e2505A ;
2013PhRvD..88i3004A  as,

$\begin{split}\nu\,(\bar{\nu})+{}^{16}\mathrm{O}&\to\nu\,(\bar{\nu})+{}^{15}\mathrm{O}^{*}+n,\\\
\nu\,(\bar{\nu})+{}^{16}\mathrm{O}&\to\nu\,(\bar{\nu})+{}^{15}\mathrm{N}^{*}+p,\end{split}$
(1)

where an excited nucleus is left in some cases that will relax with an
emission of $\gamma$-ray(s), e.g., a branching fraction of $\sim$65% is
predicted for an excited state after the NCQE interaction in Ref.
2013PhRvD..88i3004A . Energies of the knocked-out nucleon are typically in
hundreds of MeV. Such nucleons, particularly neutrons, would cause additional
reactions inside the detector that produce more $\gamma$-rays, as
schematically shown in Figure 1. Here, the primary $\gamma$-ray from the NCQE
interaction and the secondary $\gamma$-ray from the following nuclear
reactions are emitted within only O(10) ns, which makes it difficult to
separate from each other, and then both $\gamma$-rays are treated as signals
in analysis. For this reason, an accurate modeling of secondary $\gamma$-ray
emission in the neutron-oxygen reaction is essential to predict the signal
distribution in the NCQE measurement at water Cherenkov detectors. Currently,
it is implied that there is an issue with modeling of the secondary
$\gamma$-ray emission, causing a large systematic uncertainty in the cross
section measurement 2019PhRvD.100k2009A ; 2024PhRvD.109a1101S . In the T2K
measurement 2019PhRvD.100k2009A , more events are expected at large Cherenkov
angles, where multiple $\gamma$-rays are typically reconstructed, in the Monte
Carlo simulation than the observed data. In the SK measurement
2024PhRvD.109a1101S , several nuclear models were tested for the NCQE sample
and a potential ill-modeling of the secondary $\gamma$-ray emission is pointed
out. Therefore, there is a high demand that more experimental data are
provided for better modeling of the nuclear $\gamma$-ray production in water.

Figure 1: Schematic illustration of the neutrino-oxygen NCQE interaction and
the secondary interaction by knocked-out neutrons.

The ${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross sections have been measured
for various $\gamma$-ray components using quasi-mono energetic neutron beams
at Osaka University’s Research Center for Nuclear Physics (RCNP). In order to
cover the energy range of the knocked-out nucleon induced by the neutrino-
oxygen NCQE interaction, neutron energies of 30, 80, and 250 MeV were chosen
for the measurements. The first RCNP experiment was conducted in 2017 as E487
at $E_{n}=80$ MeV 2024PhRvC.109a4620A . Additional beam experiments, E525,
were performed at $E_{n}=30$ and 250 MeV, and their
${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross section results are reported in
this paper. The rest part of the paper is structured as follows. In Section 2,
the RCNP facility and our experimental setup are described. In Sections 3 and
4, the neutron flux and $\gamma$-ray measurements are detailed. The
${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross section results from the
experiments are reported in Section 5. Discussion about our results and
impacts to the NCQE measurement at water Cherenkov detectors is given in
Section 6, before concluding in Section 7.

##  2 The RCNP-E525 Experiment

The E525 experiment was carried out in a 100 m long beamline at RCNP. The
basic setup and measurement strategy followed the previous experiment (E487)
performed in the same beamline 2024PhRvC.109a4620A , while several parts have
been optimized or improved. Protons were accelerated by two cyclotrons, the
K140 AVF cyclotron and the K400 ring cyclotron, and guided into a lithium
target to produce a quasi-mono energetic neutron beam via the ${}^{7}{\rm
Li}(p,n){}^{7}{\rm Be}$ reaction. The lithium target consists of 7.5%
${}^{6}{\rm Li}$ and 92.5% ${}^{7}{\rm Li}$ and its thickness is 10.0$\pm$0.2
mm. In E525, we performed measurements at two proton energy settings, 30 MeV
and 250 MeV, on December 16–17 and October 30–31 in 2018, respectively. The
resulting beam entering the experimental hall is dominated by neutral
particles as charged particles are rejected after the lithium target by a
magnetic field. The proton beam current was being monitored by a Faraday cup
placed at the beam dump and its data was read out by a scaler module.

Figure 2: Schematic drawing of the facility and beamline for the E525
experiment at RCNP. Our detectors were placed at 11 (20) m downstream of the
lithium target in the 30 (250) MeV experiment.

Figure 2 shows a schematic view of the beamline as well as our detector setup
therein. In the 30 MeV experiment, we placed an acrylic container which
measures 20 cm in diameter and 26.5 cm in length on the beam axis 11 m away
from the lithium target. In the 250 MeV experiment, in order to measure the
time-of-flight distribution accurately, we set up the water sample and
detectors farther, 20 m away from the target. Accordingly, to gain enough data
statistics, we used two identical acrylic containers and placed them as their
side facing the beam axis (see Figure 2). A high-purity germanium detector
(HPGe) and a CsI(Tl) scintillator were placed at upstream of the acrylic
container(s) off the beam axis in order to avoid neutron background. The
location was determined by looking at the neutron beam profile measured in
each experiment as described in Section 4. A liquid scintillator (BC-501A) was
aligned at downstream of the container and on the beam axis. These
container(s) and detectors are the same ones as the E487 experiment
2024PhRvC.109a4620A . All detector waveforms were recorded with a 14-bit
250-MHz CAEN DT5725 digitizer. Measurements were performed in two
configurations with the container(s) being water-filled and empty.

There are several updated aspects in E525 from the previous E487 experiment.
The main $\gamma$-ray detector in E525 is HPGe whose energy resolution is much
better than the ${\rm LaBr_{3}(Ce)}$ scintillator used in E487. This helped us
to distinguish multiple $\gamma$-ray peaks more clearly. In E487, there might
be a potential issue with the proton current measurement by a Faraday cup and
a scaler module. The proton current setting was changed frequently back-and-
forth between the neutron flux measurement ($\sim$1 nA) and the $\gamma$-ray
measurement ($\gtrsim$100 nA) because of the setup and this required a
frequent switch of the scaler’s scale setting, which might lead to a potential
human error causing a wrong normalization in the $\gamma$-ray measurement.
This will be described in more detail in Section 6. We took a special care of
this issue in the 250 MeV experiment by avoiding frequent switch of the beam
power. In the 30 MeV experiment, in addition, we monitored the neutron flux
during the $\gamma$-ray measurement using the BC-501A detector placed at
downstream of the water sample.

##  3 Neutron Flux

Neutron flux measurement was performed using a liquid scintillator, BC-501A,
and repeated between the main $\gamma$-ray measurements in order to monitor
the beam stability throughout the experiment. The water sample was removed
during the flux measurement in order to avoid attenuation by water. In the
following part, we describe each essential portion in the flux analysis: a
pulse shape discrimination (PSD) technique for neutron identification, a time-
of-flight (ToF) method for reconstruction of the neutron kinetic energy, and
the SCINFUL-QMD Monte Carlo code bib:SCINFUL_JAEA  used for neutron detection
efficiency by BC-501A.

###  3.1 Neutron event selection

Scintillation light emission in an organic scintillator from nuclear recoil,
e.g., a proton recoiled by neutron elastic scattering, contains a larger
component with a slow time constant than electron recoil due to larger
$dE/dx$. This phenomenon appears in the shape of the resulting waveform, i.e.,
neutron-like events show a longer tail than $\gamma$-ray-induced events. In
order to characterize the waveform shape, in particular the fraction of the
tail part, the PSD parameter is defined as,

$PSD=\frac{Q_{\mathrm{tail}}}{Q_{\mathrm{total}}},$ (2)

where $Q_{\mathrm{total}}$ and $Q_{\mathrm{tail}}$ are the integrated charge
for the entire waveform and its tail part, respectively. We optimized the
integral time windows for $Q_{\mathrm{total}}$ and $Q_{\mathrm{tail}}$ as
$\left[T_{0},T_{0}+600~{}\mathrm{ns}\right]$ and
$\left[T_{0}+48,T_{0}+600~{}\mathrm{ns}\right]$, respectively, where $T_{0}$
is defined as a rise time of the waveform. The left panel of Figure 3 shows a
PSD parameter as a function of $Q_{\mathrm{total}}$ from the 30 MeV
experiment. Here, our optimized criteria for selecting neutron-like events,
$PSD>0.19$, is also shown. Because the saturation effect of electronics gets
more visible at high charges and was actually observed in the 250 MeV
experiment, the selection curve is modified as,

$PSD>p_{0}+p_{1}\times Q_{\mathrm{total}}+p_{2}\times Q_{\mathrm{total}}^{2}.$
(3)

With optimization, we determined $p_{0}=0.181$, $p_{1}=-5.79\times 10^{-4}$,
and $p_{2}=1.04\times 10^{-5}$, as shown in the right panel of Figure 3. The
selection efficiencies of the PSD cut achieved more than 99% in both
experiments, which gives a negligibly small effect on neutron detection
efficiency discussed later. Mis-identification ratio of electron recoil event
to nuclear recoil in the criterion was estimated as 0.8% in the 30 MeV
experiment and 0.5% in the 250 MeV experiment.

Figure 3:  Calculated PSD parameters as a function of $Q_{\mathrm{total}}$ for
the 30 MeV (left) and 250 MeV experiment (right). The cut criteria for neutron
selection are displayed in red. Due to the saturation effect of electronics, a
correlation between $Q_{\mathrm{total}}$ and
$PSD=Q_{\mathrm{tail}}/Q_{\mathrm{total}}$ shows a non-linearity in the large
charge region in the 250 MeV experiment.

###  3.2 Kinetic energy reconstruction

Neutron kinetic energies are reconstructed based on ToF, a time interval
between the rise time of the BC-501A waveform and a RF-synchronized periodic
pulse timing. Figure 4 shows the ToF distributions obtained in each
experiment. The sharp peaks around 120 ns in the left panel and 370 ns in the
right panel correspond to events of $\gamma$-rays instantly emitted due to the
proton-lithium interaction (called “flash $\gamma$-rays”). One can find
another peak $\sim$7 ns before the 370 ns peak in the 250 MeV experiment. From
the detailed investigation, we confirmed that a polyimide tape inserted in the
proton beam pipe for preparation of a preceding experiment caused this peak.
The position where the polyimide tape was inserted was approximately 2.5 m
upstream of the lithium target, which is consistent with the observed time
difference of $\sim$7 ns. In addition, there were no material like the
polyimide tape at the time when we performed the 30 MeV experiment. Thus, we
identified the latter peak around 370 ns as the flash $\gamma$-ray peak of the
proton-lithium interaction in the 250 MeV experiment. The resulting kinetic
energy distributions converted from ToF are shown in Figure 5. The peak
energies, 18–30 MeV and 235–270 MeV, reconstructed in each experiment are
consistent with the proton beam settings.

Figure 4: ToF distributions categorized as $\gamma$-like (red), neutron-like
(green) and all events (black) in the 30 MeV (left) and 250 MeV experiment
(right). The peak at early times mainly consisting of the $\gamma$-like events
corresponds to the flash $\gamma$-ray. The neutron from the ${}^{7}{\rm
Li}(p,n){}^{7}{\rm Be}$ reaction forms the peak at late times. Figure 5:
Reconstructed neutron kinetic energy distributions based on ToF in the 30 MeV
(left) and 250 MeV experiment (right).

###  3.3 Detection efficiency

A dedicated Monte Carlo code for computing neutron detection efficiency of the
BC-501A scintillator, SCINFUL-QMD bib:SCINFUL_JAEA , is used for the
efficiency calculation. It deals with the process of neutron reaction in
organic scintillation material and conversion to PMT light yields. Inputs to
SCINFUL-QMD include detector and source geometries, detector threshold,
attenuation factor of the scintillation material and the PMT response
function. The detector threshold is confirmed using the Compton edge of the
4.4 MeV $\gamma$-ray from an ${\rm{}^{241}Am/Be}$ source. The scintillator
attenuation factor is set to 0.008 $\rm{cm}^{-1}$ following Ref.
2024PhRvC.109a4620A . Three functional forms are available as the PMT light
output model and the difference of the resulting detection efficiency from
these three models is taken into account as a systematic uncertainty in this
analysis. Figure 6 shows the calculated neutron detection efficiency of
BC-501A as a function of neutron energy with a linear interpolation for the
energy range of 80–150 MeV. Note that because SCINFUL-QMD unnatually performs
its calculation using a constant cross section value in 80 to 150 MeV, its
developer recommends the interpolation prescription in this range.

Figure 6: Neutron detection efficiencies of the BC-501A scintillator. Red
points correspond to the SCINFUL-QMD calculations and a black line gives an
interpolation between 80 and 150 MeV as recommended by the developer.

###  3.4 Flux estimation

Neutron flux distributions in each experiment are obtained by scaling the
reconstructed neutron kinetic energy distributions with the detection
efficiency and the solid angle as shown in Figure 7. The estimated neutron
flux in the peak region is $2.11\times 10^{10}~{}\rm{neutrons/sr/\mu C}$
(18–30 MeV) and $5.01\times 10^{9}~{}\rm{neutrons/sr/\mu C}$ (235–270 MeV),
respectively. The flux results will be used when calculating the
${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross sections. Table 1 summarizes
uncertainties in the estimation of neutron flux within the 18-30 MeV and
235-270 MeV energy ranges. The statistical uncertainties are less than 1%.
Each source of systematic uncertainty is described below.

Figure 7: Neutron flux distributions in the unit of injected proton beam
current and detector solid angle in the 30 MeV experiment (left) and the 250
MeV experiment (right).

#### Detection efficiency

Ref. bib:BC-501A_attenuation  reports an uncertainty derived from physics
model used in SCINFUL-QMD is 10.0% and 15.0% for the neutron energy below and
above 80 MeV, respectively. The effect of the SCINFUL inputs, e.g., a
threshold value, light attenuation factor, choice of the light output
function, and Monte Carlo statistical uncertainty, are all smaller than 0.1%.
Consequently, 10.0% and 15.0% uncertainties are assinged to the detection
efficiency for the 30 MeV experiment and the 250 MeV experiment, respectively.

#### Beam stability

As mentioned above, the neutron flux measurements were performed between
$\gamma$-ray measurements for monitoring stability. It was stable within 1.1%
and 10.0% for the 30 MeV and 250 MeV experiment, respectively. They are
considered as the systematic uncertainty for the neutron flux. It is natural
that there should be a larger fluctuation in the 250 MeV experiment because an
electromagnet of the ring cyclotron did not reach at thermal equilibrium
sufficiently. Because no other experiments were scheduled before our beam
time, beam operation for the 250 MeV experiment began with cold magnet. On the
other hand, the AVF cyclotron was in operation before the 30 MeV experiment
and already reached at thermal equilibrium at that time; thus, the beam was
stable in the 30 MeV experiment.

#### Neutron selection with PSD

The systematic uncertainty associated with the PSD cut is estimated based on
contamination of the $\gamma$-like events in the neutron-like events. We
consider 0.8% and 0.5% as an uncertainty for each 30 and 250 MeV experiment,
respectively.

#### Solid angle

The uncertainty of the geometrical setup is estimated by conservatively
assuming the scale was misread by 1 cm. The result is 0.2% in the 30 MeV
experiment and 0.1% in the 250 MeV experiment.

Table 1: Statistical and systematic uncertainties of the neutron flux measurement. Error Source | Size (18–30 MeV) [%] | Size (235–270 MeV) [%]
---|---|---
Statistical | 1.0 | 0.4
Detection efficiency by SCINFUL-QMD | 10.0 | 15.0
Beam stability | 1.1 | 10.0
Neutron selection | 0.8 | 0.5
Solid angle | 0.2 | 0.1
Total | 10.1 | 18.0

##  4 Nuclear Deexcitation $\gamma$-rays

###  4.1 Detector positioning

Figure 8 shows the neutron-like event rates measured using BC-501A at some
distances horizontally away from the beam center in the 30 MeV and 250 MeV
experiments. The measurement method follows the one outlined in Section 3.
From this measurement, we decided to place the HPGe detector at 24 cm away
from the beam center, where the neutron event rate is reduced by an order of
$\sim$3 and $\sim$1 compared to the center in the 30 MeV and 250 MeV
experiment, respectively. A wider beam profile in the 250 MeV experiment is
because our setup is made at more downstream as explained in Section 2. We
also placed the CsI(Tl) scintillator at the same distance from the beam center
as the HPGe detector.

Figure 8: Neutron beam profile at each distance from the beam center in the 30
MeV (left) and 250 MeV experiment (right). The event counting is performed for
the peak region, 18–30 MeV for the 30 MeV experiment and 235–270 MeV for the
250 MeV experiment.

###  4.2 The HPGe waveform analysis

We designed an HPGe waveform analysis to reconstruct the time and energy for
each detected event. Figure 9 shows an example HPGe waveform taken by the
digitizer. The baseline is calculated for each waveform as an average of the
ADC counts for the range of 0 to 300 channels, where 1 channel corresponds to
4 ns in time. Here, the standard deviation ($\sigma$) is computed together and
used for determining the detected time. An event time is reconstructed as an
intersection of the calculated baseline and the waveform rise edge which is
obtained from a linear fitting to the points of 5$\sigma$ to 12$\sigma$ away
from the baseline, as schematically shown in Figure 9. For the energy
reconstruction, we use the maximum height from the baseline in each waveform,
which is justified by the fact that the pre-amplifier output from HPGe is
proportional to the waveform size.

Figure 9: An example waveform from the HPGe detector (left) and a schematic
drawing of the analysis scheme (right). In the left panel, 1 channel
corresponds to 4 ns.

The defined physics variables were calibrated before each 30 and 250 MeV
experiment on its energy and timing measurements. For the energy calibration,
we made use of $\gamma$-rays from various radioactive isotopes, including
${}^{2}\rm{H}$, ${}^{60}\rm{Co}$, ${}^{40}\rm{K}$, ${}^{241}\rm{Am/Be}$, and
${}^{56}\rm{Fe}$, which enabled us to calibrate the detector up to $\sim$8 MeV
with making sure its absolute scale within a 1% level. The energy scale was
monitored by performing a quick calibration with a few sources during the
experiments and no significant shift was observed. We also calibrated the time
response of the HPGe detector within a few ns uncertainty using two
$\gamma$-rays emitted at the same time from a ${}^{60}\rm{Co}$ source with
tagging one using a fast response detector, ${\rm LaBr_{3}(Ce)}$, and the
other by HPGe

###  4.3 Analysis time window

We performed the ToF analysis for the HPGe detector similarly to the BC-501A
data as described in Section 3. The resulting ToF distributions are shown in
Figure 10. In the distribution from the 30 MeV experiment, a clear peak by the
flash $\gamma$-ray is observed around 250 ns; however, it is not clearly
observed in the 250 MeV experiment. This is perhaps because the background
amount is expected to be more as is implied by the wider neutron beam profile
in Figure 8 and the higher energy beam. Therefore, we decide to utilize the
BC-501A data to find the flash $\gamma$-ray point which is the baseline point
for kinetic energy reconstruction. For this, we smeared the BC-501A’s ToF
distribution by resolution (standard deviation) difference in the time
reconstruction between HPGe (6.8 ns) and BC-501A (1.9 ns). The smeared ToF
distribution of the BC-501A detector is also shown in the right panel of
Figure 10.

We define two time windows, the on- and off-timing window, based on the
reconstructed ToF distributions, which are used later for obtaining the signal
and background energy spectra. Each time window is indicated in Figure 10. The
on-timing window corresponds to the range of 18–30 MeV and 235–270 MeV in
neutron energies in the 30 and 250 MeV experiment, respectively. The off-
timing window is defined using the ToF range where neutron rate is expected to
be low. It is the range between the flash $\gamma$-ray peak and the bulk
neutrons in the 30 MeV experiment, and defined as the range preceding the
expected flash $\gamma$-ray position in the 250 MeV experiment, as indicated
in Figure 10.

Figure 10: ToF distributions of the HPGe detector (blue) in the 30 MeV (left)
and 250 MeV experiment (right). In the right panel, the smeared ToF
distribution of the BC-501A detector (red) is also shown. The on-timing and
off-timing windows for each experiment are displayed as the arrows in magenta
and orange, respectively.

###  4.4 Observed $\gamma$-ray peaks

Figure 11 shows the measured energy spectra in two configurations and for the
on- and off-timing windows from the 30 MeV and 250 MeV experiments. Note that
the spectra for the off-timing window is normalized to the size of the on-
timing window. Multiple $\gamma$-ray peaks are observed in the water-filled
configuration. The observed peaks, their parent nuclei and likely physics
processes are summarized in Table 2. Many of them were observed also in the
previous experiment (E487) at $E_{n}=80$ MeV 2024PhRvC.109a4620A . The
$\gamma$-rays of 7.12, 6.92, 6.13, and 2.74 MeV are expected from different
excited states of ${\rm{}^{16}O}$ and clearly observed in both experiments.
The parent process for these $\gamma$-rays, the $(n,n^{\prime})$ scattering,
does not cause an additional nucleon in the final state. The 6.32, 5.27, and
2.30 MeV $\gamma$-rays are emitted from excited states of ${\rm{}^{15}N}$,
where a proton is released from the $(n,np)$ scattering. Here, the final state
neutron and proton may form a bound state of a deuteron, but our measurement
is not able to tell this. It should also be note that the proton release
possibly occurs after the $(n,n^{\prime})$ scattering via
${\rm{}^{16}O^{*}}\to{\rm{}^{15}N^{*}}+p$ and again the current experiment is
not sensitive enough to distinguish this. Similarly, the 6.18 and 5.18 MeV
$\gamma$-rays are produced from excited states of ${\rm{}^{15}O}$ with an
additional neutron emission. The 5.10, 4.91, and 2.31 MeV $\gamma$-rays are
expected from excited states of ${\rm{}^{14}N}$. The 5.10 and 4.91 MeV peaks
are not observed very clearly, but they are considered in the later spectral
fitting. The 2.31 MeV peak cannot be distinguished from the 2.30 MeV
$\gamma$-ray peak by the HPGe detector, hence we measure the inclusive cross
section for these two $\gamma$-rays. Intense spectral peaks are formed by the
3.68 MeV $\gamma$-ray from ${\rm{}^{13}C}(\frac{3}{2}^{-})$ and the 4.44 MeV
$\gamma$-ray from ${\rm{}^{12}C}(2^{+})$. Many peaks originating from thermal
neutron capture process are also observed: 2.22 MeV from ${\rm{}^{2}H}$, 3.84
MeV from ${\rm{}^{17}O}$, and 7.63 MeV from ${\rm{}^{57}Fe}$. In addition, the
1.46 MeV $\gamma$-ray from ${\rm{}^{40}K}$ and the 2.61 MeV $\gamma$-ray from
${\rm{}^{208}Tl}$ are produced in the highly-radiated environment. These
$\gamma$-rays are not the measurement targets of the present work.

Figure 11: Observed energy distributions in two configurations on the acrylic container(s) and for two timing windows from the 30 MeV (top) and 250 MeV experiment (bottom). Table 2: Summary of the neutron-oxygen reaction-induced $\gamma$-rays identified by the HPGe detector with their parent excited nuclei and physics processes. Energy [MeV] | Parent nucleus ($J^{\pi}$) | Physics process
---|---|---
2.30 | ${\rm{}^{15}N}(\frac{7}{2}^{+})$ | ${\rm{}^{16}O}(n,np){\rm{}^{15}N^{*}}$
2.31 | ${\rm{}^{14}N}(0^{+})$ | ${\rm{}^{16}O}(n,2np){\rm{}^{14}N^{*}}$
2.74 | ${\rm{}^{16}O}(2^{-})$ | ${\rm{}^{16}O}(n,n^{\prime}){\rm{}^{16}O^{*}}$
3.68 | ${\rm{}^{13}C}(\frac{3}{2}^{-})$ | ${\rm{}^{16}O}(n,\alpha){\rm{}^{13}C^{*}}$
4.44 | ${\rm{}^{12}C}(2^{+})$ | ${\rm{}^{16}O}(n,n\alpha){\rm{}^{12}C^{*}}$
4.91 | ${\rm{}^{14}N}(0^{-})$ | ${\rm{}^{16}O}(n,2np){\rm{}^{14}N^{*}}$
5.10 | ${\rm{}^{14}N}(2^{-})$ | ${\rm{}^{16}O}(n,2np){\rm{}^{14}N^{*}}$
5.18 | ${\rm{}^{15}O}(\frac{1}{2}^{+})$ | ${\rm{}^{16}O}(n,2n){\rm{}^{15}O^{*}}$
5.27 | ${\rm{}^{15}N}(\frac{5}{2}^{+})$ | ${\rm{}^{16}O}(n,np){\rm{}^{15}N^{*}}$
6.13 | ${\rm{}^{16}O}(3^{-})$ | ${\rm{}^{16}O}(n,n^{\prime}){\rm{}^{16}O^{*}}$
6.18 | ${\rm{}^{15}O}(\frac{3}{2}^{-})$ | ${\rm{}^{16}O}(n,2n){\rm{}^{15}O^{*}}$
6.32 | ${\rm{}^{15}N}(\frac{3}{2}^{-})$ | ${\rm{}^{16}O}(n,np){\rm{}^{15}N^{*}}$
6.92 | ${\rm{}^{16}O}(2^{+})$ | ${\rm{}^{16}O}(n,n^{\prime}){\rm{}^{16}O^{*}}$
7.12 | ${\rm{}^{16}O}(1^{-})$ | ${\rm{}^{16}O}(n,n^{\prime}){\rm{}^{16}O^{*}}$

###  4.5 Background estimation

For the $\gamma$-ray production cross section measurement to be performed
later, we estimate background events to extract the intensity of $\gamma$
signals. In this study, four types of background events are considered as
following the previous experiment 2024PhRvC.109a4620A : fast neutron hits,
non-water background, time-uniformly distributed $\gamma$ and $\beta$
background, and scattered neutron-induced events. Each of them is described in
the following part.

#### Fast neutrons directly entering the detector

Neutrons which are scattered by the acrylic container(s) or incident off beam
axis and enter the HPGe directly can be background events. For the measurement
of this background, a CsI(Tl) scintillator was placed at upstream of the
acrylic container in a symmetric position with respect to HPGe across the beam
axis. Based on the PSD analysis on CsI(Tl) data, which has been established in
the previous tests 2018PTEP.2018d3H01A ; 2024PhRvC.109a4620A , such background
is negligible up to $\sim$7 MeV and accounts for 2% of the observed
$\gamma$-ray intensity in the 7–8 MeV region.

#### Non-water background

$\gamma$-rays produced from sources other than water, such as the acrylic
container(s) and experimental stage, constitute background. This can be
subtracted by the measurement with an empty configuration, as shown as the
blue histogram in Figure 11.

####  $\gamma$-rays from thermal neutron capture and $\beta$’s from beta decay

$\gamma$-rays originating from thermal neutron capture and $\beta$’s from beta
decay are expected to distribute uniformly in time as their emission timescale
is much longer than the beam repetition cycle. Therefore, we estimate these
background events by the off-timing spectrum created with a TOF cut (see
Figure 11).

#### Scattered fast neutron-induced $\gamma$-rays

There still remains continuous background in the observed HPGe spectrum, which
are expected to arise from reactions between neutrons scattered off water and
the surrounding materials. Such background was confirmed in E487
2024PhRvC.109a4620A  and its estimation method was established. Therefore, we
take the same estimation method in the current study, which will be described
as a continuous background template in Section 5.

##  5 Cross Section Analysis

In this section, we describe the ${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross
section measurement for the peaks shown in Table 2. We take a binned fitting
method on the energy spectrum obtained after some background subtraction.
First, we subtract the non-water background and the uniform-in-time background
(see Section 4.5) as,

$N^{\mathrm{Data}}_{\mathrm{bin}}=\left(N^{\mathrm{Water,ON}}_{\mathrm{bin}}-N^{\mathrm{Water,OFF}}_{\mathrm{bin}}\right)-\left(N^{\mathrm{Empty,ON}}_{\mathrm{bin}}-N^{\mathrm{Empty,OFF}}_{\mathrm{bin}}\right),$
(4)

where $N^{\mathrm{Water,ON}}_{\mathrm{bin}}$ is the number of events in the
bin-th energy bin of the observed spectrum from the water-filled container
data in the on-timing window. Accordingly,
$N^{\mathrm{Water,OFF}}_{\mathrm{bin}}$ corresponds to that from the water-
filled container data in the off-timing window.
$N^{\mathrm{Empty,ON}}_{\mathrm{bin}}$ and
$N^{\mathrm{Empty,OFF}}_{\mathrm{bin}}$ correspond to the empty-container data
for the on- and off-timing window, respectively. We then perform the fitting
to this spectrum after these subtraction with including the remaining
continuous background as a part.

###  5.1 Signal and background templates

Spectrum templates of $\gamma$-rays from discrete levels of excited nuclei are
produced by a Monte Carlo simulation based on the GEANT4 package bib:geant4
in order to include a geometrical effect of the setup on the spectra observed
in the HPGe detector. The simulated $\gamma$-rays are listed in Table 2. The
$\gamma$-rays with isotropic momentum direction are generated at a point
following an emission point probability density function (PDF) in the water
sample. The emission point PDF was generated based on the position where
neutron inelastic scattering occurred in the MC simulation shooting neutron
with the measured beam profile (Fig. 8) onto the water target. There is almost
no change in the spectrum template if the PDF is replaced with a uniform
distribution. Thus, a systematic uncertainty due to the emission point is
negligibly small.

The Doppler broadening effect is taken into consideration for the states with
a shorter lifetime than their interaction time scale. Thus, we apply a 1%
smearing with a Gaussian distribution for the 2.30, 2.74, 3.68, 4.91, 5.18,
6.18, 6.32, 6.92 and 7.12 MeV $\gamma$-rays, and gives a 2% smearing for the
4.44 MeV $\gamma$-ray, respectively. Because of the Doppler broadening and the
detector resolution, it is difficult to separate contributions from the 2.30
MeV and 2.31 MeV $\gamma$-rays. Thus, the 2.30 MeV $\gamma$-ray template is
included in the fit for the inclusive estimate of these two peaks.

The simulated deposited energy in the HPGe detector is smeared with an energy
resolution function. It is determined by fitting energy resolution at
calibration points with a function $\sigma/E=p_{0}+p_{1}/E$
($\sigma/E=p_{0}+p_{1}/E^{2}$) for the 30 (250) MeV experiment, where $\sigma$
is the standard deviation determined by a Gaussian fit and $E$ is the peak
energy. As a result of the curve fitting, the parameters, ($p_{0},p_{1}$), are
determined as ($(3.51\pm 0.03)\times 10^{-1},(201.5\pm 2.8)\times 10^{3}$) for
the 30 MeV experiment and ($(2.32\pm 0.10)\times 10^{-1},(5.38\pm 0.20)\times
10^{2}$) for the 250 MeV experiment, respectively. An effect of the fitting
errors on the cross section estimation is negligible.

An uncertainty about the $\gamma$-ray detection efficiency at the HPGe
detector is estimated based on the calibration data in which a 60Co source was
placed at 6, 9 and 12 cm away from the HPGe detector. Simulated spectra in the
identical setup with the calibration were compared with the calibration data.
We found that the simulated spectra reproduced the observed data within a 10%
difference at each source position. Thus, this 10% is taken into consideration
as a systematic uncertainty in the spectrum fitting.

As discussed above, the continuous background spectrum formed by $\gamma$-rays
from the scattered neutron reactions on surrounding materials is considered in
the fit. Its spectrum shape is estimated based on the on-timing data with the
empty acrylic container(s). The template as a smooth function is obtained by
fitting the empty sample data with an exponential function
$f(E)=\exp(p_{0}+p_{1}E)$, where $E$ is the observed energy at the HPGe
detector in the unit of keV. The determined parameters are $p_{0}=8.67\pm
0.06$ and $p_{1}=(-5.52\pm 0.15)\times 10^{-4}$ for the 30 MeV experiment, and
$p_{0}=9.65\pm 0.05$ and $p_{1}=(-5.27\pm 0.12)\times 10^{-4}$ for the 250 MeV
experiment. The obtained spectra from two experiments are consistent to each
other within 10%.

###  5.2 Least chi squared fit

The data spectrum is fit with the simulated templates based on a chi-squared
taking the systematic uncertainty into consideration. The definition of the
chi-squared is denoted as,

$\begin{split}\chi^{2}&=\sum_{\mathrm{bin}}\left(\frac{N^{\mathrm{Data}}_{\mathrm{bin}}-\eta\times
N^{\mathrm{MC}}_{\mathrm{bin}}}{\sigma^{\mathrm{stat}}_{\mathrm{bin}}}\right)^{2}+\left(\frac{1-\eta}{\sigma^{\mathrm{syst}}}\right)^{2},\\\
&N^{\mathrm{MC}}_{\mathrm{bin}}=\sum_{j}f_{j}\times
n_{\mathrm{bin}}^{j},\end{split}$ (5)

where $\eta$ is a nuisance parameter characterizing the systematic uncertainty
from the $\gamma$-ray detection efficiency at the HPGe detector. As mentioned
above, we assigned a 10% uncertainty on the detection efficiency.
$N^{\mathrm{Data}}_{\mathrm{bin}}$ and $N^{\mathrm{MC}}_{\mathrm{bin}}$
represent the number of the observed events at each energy bin for the
observed data and the Monte Carlo simulation, respectively. The Monte Carlo
prediction, $N^{\mathrm{MC}}_{\mathrm{bin}}$, consists of the summation of the
spectrum templates ($n^{j}_{\mathrm{bin}}$) with scale factors ($f_{j}$) where
$j$ corresponds to the template number, 0: continuous background, 1: 2.30 &
2.31 MeV, 2: 2.74 MeV, 3: 3.68 MeV, 4: 4.44 MeV, 5: 4.91 MeV, 6: 5.10 MeV, 7:
5.18 MeV, 8: 5.27 MeV, 9: 6.13 MeV, 10: 6.18 MeV, 11: 6.32 MeV, 12: 6.92 MeV,
and 13: 7.12 MeV. The statistical errors of the observed data and simulation
are considered in $\sigma^{\mathrm{stat}}_{\mathrm{bin}}$ at each energy bin.
Minimization of the chi-squared and error estimation were performed using a
dedicated tool, MINUIT James:1975dr , with the scale factors and the nuisance
parameter as fitting parameters. The best-fit spectra are shown in Figures 12
and 13, and the obtained scale factors are summarized in Table 3. We show the
best-fit results with $\pm 1\sigma$, but only give the 90% upper limits for
some components whose scale factor is zero-consistent as a result of the fit.
The minimum chi-squared value is 564.2 (135.7) in the degree of freedom 635
(89) for the 30 (250) MeV case. Therefore, the assumption that the observed
spectra consist of $\gamma$-rays from discrete levels of nuclei and the
continuum component works well in terms of the goodness of fit.

Figure 12: Energy spectra of the signal $\gamma$-ray in the 30 MeV experiment
(black points) overlaid with the best-fit of the spectrum fitting (dashed red
line) in a linear scale (top) and a logarithmic scale (bottom). The signal and
background templates are overlaid to show the result of fit, which are colored
in as follows: orange (7.12 MeV), azure (6.92 MeV), green (6.32 MeV), magenta
(6.18 MeV), light blue (6.13 MeV), pink (5.27 MeV), blue (5.18 MeV), light
green (4.44 MeV), light orange (3.68 MeV), yellow (2.74 MeV), purple
(2.30/2.31 MeV) and gray (background continuum).

Figure 13: Energy spectra of the signal $\gamma$-ray in the 250 MeV experiment
(black points) overlaid with the best-fit of the spectrum fitting (dashed red
line) in a linear scale (top) and a logarithmic scale (bottom). The color
notation is the same as Figure 12.

###  5.3 $\gamma$-ray production cross section

The $\gamma$-ray production cross section at each discrete level can be
computed as,

$\sigma_{j}^{\gamma}=f_{j}\times\frac{N^{\mathrm{MC}}_{\mathrm{gen}}}{T_{\mathrm{Oxygen}}\times\phi_{n}},$
(6)

where $N^{\mathrm{MC}}_{\mathrm{gen}}$ is the number of generated
$\gamma$-rays in producing the spectrum template, $T_{\mathrm{Oxygen}}$ is the
number of 16O existing in the water sample, and $\phi_{n}$ is the neutron
flux. The numbers used here are as follows:
$N^{\mathrm{MC}}_{\mathrm{gen}}=1.0\times 10^{8}$ ($1.0\times 10^{8}$),
$T_{\mathrm{Oxygen}}=2.3\times 10^{26}$ ($4.6\times 10^{26}$),
$\phi_{n}=1.4\times 10^{4}$ ($1.2\times 10^{3}$) neutrons/cm2/$\mu$C for the
30 (250) MeV experiment. Table 3 shows the summary of the measured
$\gamma$-ray production cross sections in the E525 experiment. In addition to
the error from the fit, the neutron flux uncertainty is propagated for the
total uncertainty on the cross section result (see Table 1). As mentioned
above, we give the 90% upper limit for the peaks whose statistics was not
significant to be measured in the fit. The total $\gamma$-ray production cross
sections are measured as $\sigma_{\mathrm{total}}=177.6\pm 12.3$ at 30 MeV and
$\sigma_{\mathrm{total}}=75.9\pm 9.7$ at 250 MeV by integrating over the 2 to
8 MeV $\gamma$-rays and subtracting the continuum component. Discussion about
the results is given in the next section.

Table 3: The $\gamma$-ray production cross section results from the 30 MeV and 250 MeV experiments. For the 2.30 and 2.31 MeV $\gamma$-rays, the inclusive measurement results are shown. No. | Peak Energy [MeV] | Cross Section (30 MeV) [mb] | Cross Section (250 MeV) [mb]
---|---|---|---
1 | 2.30/2.31 | $0.8\pm 0.7$ |  $<6.9$ (90% C.L.)
2 | 2.74 | $2.2\pm 0.9$ |  $<2.8$ (90% C.L.)
3 | 3.68 | $4.1\pm 1.2$ | $3.7\pm 2.6$
4 | 4.44 | $37.5\pm 6.8$ | $10.8\pm 3.3$
5 | 4.91 | $1.9\pm 1.2$ |  $<3.7$ (90% C.L.)
6 | 5.10 |  $<0.1$ (90% C.L.) | $5.4\pm 2.6$
7 | 5.18 | $1.3\pm 0.9$ |  $<6.2$ (90% C.L.)
8 | 5.27 | $9.9\pm 1.9$ | $8.0\pm 2.7$
9 | 6.13 | $22.3\pm 4.1$ | $9.4\pm 3.1$
10 | 6.18 | $4.6\pm 1.3$ | $9.2\pm 3.4$
11 | 6.32 | $12.9\pm 2.5$ | $15.3\pm 4.2$
12 | 6.92 | $6.0\pm 1.5$ | $7.2\pm 2.6$
13 | 7.12 | $9.0\pm 1.9$ | $5.3\pm 2.3$

##  6 Discussion

Ref. bib:nelson  reported the measurement results of the
${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross sections using a white neutron
beam which has a broad kinetic energy range up to 200 MeV. This work gives the
unique results using nearly mono energy neutrons at 30 and 250 MeV, and the
cross sections at 250 MeV are the first measurement results above 200 MeV. In
addition, more types of $\gamma$-rays are identified in this work than the
work done in Ref. 2024PhRvC.109a4620A  because of the improvement on the
energy resolution thanks to the HPGe detector. For example, the $\gamma$-rays
of 2.74 MeV from $\rm{}^{16}O(2^{-})$, 4.91 MeV from $\rm{}^{14}N(0^{-})$,
5.10 MeV from $\rm{}^{14}N(2^{-})$, and 6.18 MeV from
$\rm{}^{15}N(\frac{3}{2}^{-})$ are newly measured. It is difficult to identify
the 5.10 MeV $\gamma$-ray based on its photopeak in Ref. bib:nelson  suffering
from the double escape peak of the intense 6.13 MeV $\gamma$-ray. However, our
spectrum fit method improves sensitivity to such peaks due to information of
spectrum shape, e.g., the single and double escape peaks and Compton edge.

In the present work, the 6.18 MeV and 6.32 MeV $\gamma$-rays are observed and
their cross sections are measured. They are due to deexcitation from a
$\frac{3}{2}$-hole state of $\rm{}^{15}O$ and $\rm{}^{15}N$, respectively,
which are expected to be generated as a result of the NCQE reaction by
neutrinos as well 2012PhRvL.108e2505A ; 2013PhRvD..88i3004A . The measured
cross section of the 6.32 MeV $\gamma$-ray is larger than that of the 6.18 MeV
$\gamma$-ray both at 30 and 250 MeV. This is considered because a neutron is
more likely to collide with a proton than a neutron due to the nucleon
isospin.

Figure 14: Measured cross sections of the 4.44 MeV (top left), 5.27 MeV (top
right), 6.13 MeV (bottom left) and total (bottom right) $\gamma$-ray
production as a function of neutron kinetic energy. The result of this work
(red points) and E487 (green point) are compared with Ref. bib:nelson  (black)
in each plot. The 6.13 MeV $\gamma$-ray production by proton-oxygen
interaction PhysRevC.35.1214  (blue points) is displayed as well.

Figure 14 shows comparison between measurements of the
${}^{16}\mathrm{O}(n,\mathrm{X}\gamma)$ cross sections of the 4.44 MeV, 5.27
MeV, 6.13 MeV peaks and the total $\gamma$-ray production from the present
work and Refs. bib:nelson ; 2024PhRvC.109a4620A . The measured cross section
at a 250 MeV neutron energy shows consistency with the adjacent data points
from Ref. bib:nelson  within the estimated uncertainty. In contrast, our
measurement results at a 30 MeV neutron energy show smaller cross sections by
a factor of 2 to 3 than Ref. bib:nelson . This is under investigation in more
detail using a Monte Carlo simulation, but no longer discussed here as it is
beyond the current scope.

In this work, as shown in Tables 2 and 3, the measured cross sections for the
decay mode with neutron emission ($\gamma$-rays of 2.31, 4.91, 5.10, 5.18, and
6.18 MeV) are smaller compared to others. At water Cherenkov detectors,
because protons and ions with energies of current interest are not observable,
nuclear reactions with an additional neutron emission are more crucial. Our
results may indicate that the secondary $\gamma$-ray emission is likely to be
caused by the initial neutron from the primary neutrino interaction in many
cases. In the NCQE measurement at T2K 2019PhRvD.100k2009A , as mentioned in
Section 1, the observed number of events at large Cherenkov angles, where
events with more $\gamma$-rays are reconstructed, was more than the
expectation by Monte Carlo simulation in the neutrino-dominant sample. In
contrast, this mismatch was found to be mitigated for the antineutrino-
dominant sample. The kinetic energy of knocked-out nucleons by neutrino
interactions is larger than that by antineutrinos and hence the observation in
T2K implies that the current model predicts more secondary $\gamma$ emission
for higher energy neutrons. This could be improved by the measured fact in
this work that secondary $\gamma$-rays are dominantly produced by the initial
neutron. Another useful reference is measurement of the neutron multiplicity
for neutrino interactions performed by T2K using its charged-current
(CC)-dominant sample 2019TAUP.Akutsu.T2K.Neutron ; bib:rakutsu . In this
study, the observed number of neutrons after the CC interaction was
significantly smaller than the predictions by three different neutrino
interaction models. While the results are from the CC-dominant sample, the
same secondary interaction model is used as NCQE in T2K. This study accesses
the final state neutron directly than Ref. 2019PhRvD.100k2009A , as it was
done via $\gamma$-rays, and strongly supports the indication that the current
model contributes to more secondary neutron interactions. In addition, in Ref.
2024PhRvD.109a1101S  using the atmospheric neutrino NCQE sample, the Cherenkov
angle distributions from observed data and different secondary interaction
models are compared. This disfavors the secondary interaction model based on
intranuclear cascade model by Bertini G4BERT_WRIGHT2015175  used in the two
T2K studies above. Instead, the INCL++ model INCL_PhysRevC.87.014606  is
favored in this sample, and therefore comparing the prediction from this model
and the current beam test results would be interesting.

##  7 Conclusion

The E525 experiment to measure $\gamma$-rays from the neutron-oxygen
interaction was carried out at RCNP using a quasi-mono energetic neutron beam
at 30 and 250 MeV. A high energy resolution is achieved using the HPGe
detector and helps to identify multiple $\gamma$-ray peaks. Neutron flux
measurement using the BC-501A scintillator confirms the absolute flux of
quasi-mono energitic neutrons by lithium-proton interaction and its
corresponding timing window for the signal $\gamma$-rays in the HPGe detector.
We gave a special care about the choice of proton beam current in order to
avoid a potential human error on the current measurement suspected in the
former E487 experiment. The background estimation and subtraction were done
based on the off-timing and the empty container data. A spectrum fitting
method is employed based on a chi-squared to measure the cross section for
each peak and the remaining background events simultaneously. The production
cross section of thirteen kinds of $\gamma$-rays and the total of them were
measured, and this work gives the first measurement results at a 250 MeV
neutron energy. The obtained data are useful for validation of neutron
inelastic interaction models used in the detector simulation for neutrino
event prediction. Comparing the simulation result with the present
experimental data will help improving the systematic uncertainty caused by the
neutron-oxygen interaction subsequently occurring after the neutrino NCQE
interaction in water Cherenkov detectors.

## Acknowledgments

The authors are thankful to the RCNP staff for giving beam time to the
experiment and the accelerator group for supplying stable beams. This work was
supported by Japan MEXT KAKENHI Grant Number 17J06141, 18H05537, 26400292, and
25105002.

## References

  * [1] S. Fukuda et al.  Nucl. Instrum. Meth., A501:418–462, 2003.
  * [2] K. Abe et al.  arXiv, 1109.3262 [hep-ex].
  * [3] K. Abe et al.  Diffuse supernova neutrino background search at Super-Kamiokande.  Physical Review D, 104(12):122002, 2021.
  * [4] M. Harada et al.  Search for Astrophysical Electron Antineutrinos in Super-Kamiokande with 0.01% Gadolinium-loaded Water.  Astrophys. J. Lett., 951(2):L27, 2023.
  * [5] K. Abe et al.  Measurement of the neutrino-oxygen neutral-current interaction cross section by observing nuclear deexcitation $\gamma$ rays.  Physical Review D, 90(7):072012, 2014.
  * [6] K. Abe et al.  Measurement of neutrino and antineutrino neutral-current quasielasticlike interactions on oxygen by detecting nuclear deexcitation $\gamma$ rays.  Physical Review D, 100(11):112009, 2019.
  * [7] L. Wan et al.  Measurement of the neutrino-oxygen neutral-current quasielastic cross section using atmospheric neutrinos at Super-Kamiokande.  Physical Review D, 99(3):032005, 2019.
  * [8] S. Sakai et al.  Measurement of the neutrino-oxygen neutral-current quasielastic cross section using atmospheric neutrinos in the SK-Gd experiment.  Physical Review D, 109(1):L011101, 2024.
  * [9] Artur M. Ankowski, Omar Benhar, Takaaki Mori, Ryuta Yamaguchi, and Makoto Sakuda.  Analysis of $\gamma$-Ray Production in Neutral-Current Neutrino-Oxygen Interactions at Energies above 200 MeV.  Physical Review Letters, 108(5):052505, 2012.
  * [10] Artur M. Ankowski and Omar Benhar.  Neutron knockout in neutral-current neutrino-oxygen interactions.  Physical Review D, 88(9):093004, 2013.
  * [11] Y. Ashida, H. Nagata, M. Mori, G. Collazuol, D. Fukuda, T. Horai, F. Iacob, A. Konaka, Y. Koshio, T. Nakaya, C. Nantais, T. Shima, A. Suzuki, Y. Takeuchi, H. Tanaka, R. Wendell, and T. Yano.  Measurement of $\gamma$ -ray production via the neutron-16O reaction using a 77 MeV quasimonoenergetic neutron beam.  Physical Review C, 109(1):014620, 2024.
  * [12] D. Sato et al.  JAEA-Data/Code 2006-023, 2006.
  * [13] Y. Iwamoto et al.  Nucl. Instr. and Meth. Phys. Res. A, 804, 2015.
  * [14] Y. Ashida, H. Nagata, Y. Koshio, T. Nakaya, and R. Wendell.  Separation of gamma-ray and neutron events with CsI(Tl) pulse shape analysis.  Progress of Theoretical and Experimental Physics, 2018(4):043H01, 2018.
  * [15] S. Agostinelli et al.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, 506:250 – 303, 2003.
  * [16] F. James and M. Roos.  Minuit: A System for Function Minimization and Analysis of the Parameter Errors and Correlations.  Comput. Phys. Commun., 10:343–367, 1975.
  * [17] A. Michaudon R. O. Nelson, M. B. Chadwick and P. G. Young.  High-resolution measurements and calculations of photon-production cross sections for ${}^{16}\mathrm{O}$(n,x$\gamma$) reactions induced by neutrons with energies between 4 and 200 mev.  Nuclear Science and Engineering, 138(2):105–144, 2001.
  * [18] F. L. Lang, C. W. Werntz, C. J. Crannell, J. I. Trombka, and C. C. Chang.  Cross sections for production of the 15.10-mev and other astrophysically significant gamma-ray lines through excitation and spallation of ${}^{12}\mathrm{C}$ and ${}^{16}\mathrm{O}$ with protons.  Phys. Rev. C, 35:1214–1227, Apr 1987.
  * [19] R. Akustu et al.  Study of neutrons associated with neutrino and anti-neutrino interactions on water at T2K.  Journal of Physics: Conference Series, 1468(1):012123, 2020.
  * [20] R. Akutsu.  Ph.D Thesis, The University of Tokyo, 2019.
  * [21] D.H. Wright and M.H. Kelsey.  The geant4 bertini cascade.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 804:175–188, 2015.
  * [22] A. Boudard, J. Cugnon, J.-C. David, S. Leray, and D. Mancusi.  New potentialities of the liège intranuclear cascade model for reactions induced by nucleons and light charged particles.  Phys. Rev. C, 87:014606, Jan 2013.


# Privacy-Preserving Cooperative Visible Light Positioning for Nonstationary
Environment: A Federated Learning Perspective

Tiankuo Wei<EMAIL_ADDRESS>Xiamen UniversityXiamenFujianChina and
Sicong Liu<EMAIL_ADDRESS>Xiamen UniversityXiamenFujianChina

(2022)

###### Abstract.

Visible light positioning (VLP) has drawn plenty of attention as a promising
indoor positioning technique. However, in nonstationary environments, the
performance of VLP is limited because of the highly time-varying channels. To
improve the positioning accuracy and generalization capability in
nonstationary environments, a cooperative VLP scheme based on federated
learning (FL) is proposed in this paper. Exploiting the FL framework, a global
model adaptive to environmental changes can be jointly trained by users
without sharing private data of users. Moreover, a Cooperative Visible-light
Positioning Network (CVPosNet) is proposed to accelerate the convergence rate
and improve the positioning accuracy. Simulation results show that the
proposed scheme outperforms the benchmark schemes, especially in nonstationary
environments.

Federated learning; cooperative localization; visible light positioning; data
privacy; neural network

††copyright: acmcopyright††journalyear: 2022††doi:
XXXXXXX.XXXXXXX††conference: UbiComp-ISWC ’22: Adjunct Proceedings of the 2022
ACM International Joint Conference on Pervasive and Ubiquitous Computing and
Proceedings of the 2022 ACM International Symposium on Wearable Computers;
September 11–15, 2022; Atlanta, USA††booktitle: Adjunct Proceedings of the
2022 ACM International Joint Conference on Pervasive and Ubiquitous Computing
and Proceedings of the 2022 ACM International Symposium on Wearable Computers
(UbiComp/ISWC ’22 Adjunct), September 11–15, 2022, Atlanta, USA††price:
15.00††isbn: 978-1-4503-XXXX-X/XX/XX††ccs: Information systems Location based
services††ccs: Computing methodologies Cooperation and coordination††ccs:
Human-centered computing Ubiquitous and mobile computing systems and tools

##  1\. Introduction

Nowadays, demands for indoor positioning systems are dramatically increasing
due to the ongoing explosive growth of location based services. Unlike the
outdoor environment that has been covered by the Global Navigation Satellite
System (GNSS), indoor positioning is still a challenge due to the complicated
environment (Zafari et al., 2019). Recently, with the advent of visible light
communication (VLC) and the ubiquitous deployment of light-emitting diodes
(LEDs), visible light positioning (VLP) has drawn plenty of attention from
both academia and industry (Matheus et al., 2019; Yang et al., 2016, 2017). It
has been envisioned as a promising indoor positioning technique, because of
the low hardware cost, high accuracy, unregulated spectrum, and low
electromagnetic interference (Luo et al., 2017).

Most of the research on VLP investigates the characteristics of received
optical signals, such as received signal strength (RSS), time of arrival
(TOA), and angle of arrival (AOA) (Zhuang et al., 2018; Liu et al., 2015).
These characteristics are measured and mapped to the location coordinates with
the positioning algorithms, such as triangulation and fingerprinting. Over the
past few years, learning-based methods have been introduced into indoor
positioning to further improve the positioning accuracy. For example, Van _et
al._ (Van et al., 2017) proposed an RSS-based method using weighted
$k$-Nearest Neighbors ($k$-NN), which can achieve better accuracy than
trilateration. Arfaoui _et al._ (Arfaoui et al., 2021) proposed an RSS-based
fingerprinting approach utilizing multiple layer perceptron (MLP) and
convolutional neural network (CNN) to estimate both the position and
orientation of the user equipment (UE). Lin _et al._ (Lin and Zhang, 2020)
extracted the location information from the characteristic of channel impulse
response with four structures of deep neural networks.

However, the performance of the learning-base methods may degrade in
nonstationary environments, because of the time-varying channel property
caused by obstacle movement, ambient light interference, device aging, etc. To
improve the positioning accuracy and generalization capability in
nonstationary environments, cooperative localization methods have been
investigated. For example, Keskin _et al._ (Keskin et al., 2018) designed a
novel architecture that facilitates communications among VLC receivers for
cooperative localization, and proposed a location estimation algorithm
suitable for decentralized implementation. Yang _et al._ (Yang et al., 2019a)
proposed a cooperative VLP system with the sum-product algorithm to achieve
accurate positioning in the obstacle blocking situation. However, the
position-related knowledge of the users remains to be fully extracted to
further improve the positioning performance. Moreover, privacy data of users
should be well protected, as the location data usually contains much private
information.

Therefore, in this paper, we introduce the federated learning (FL) (Konečný et
al., 2016; Yang et al., 2019b) framework and propose a privacy-preserving
cooperative VLP scheme for nonstationary environments. Specifically, we design
a cooperative VLP framework, where UEs participate in data collection and
model training. Based on this framework, we propose an FL-based VLP algorithm
aiming at jointly constructing an adaptive position estimation model without
revealing their local data to the centralized server. To further improve the
performance, we propose a Cooperative Visible-light Positioning Network
(CVPosNet) to accelerate the convergence rate and improve the positioning
accuracy.

The remainder of this paper is organized as follows. The system model of the
indoor VLP is presented and the positioning problem is formulated in Section
2. The proposed cooperative VLP scheme is introduced in Section 3. The
simulation results are reported and discussed in Section 4, followed by the
conclusions in Section 5.

##  2\. System Model and Problem Statement

###  2.1. System Model of Indoor Visible Light Positioning

Figure 1. The channel model of the indoor VLP system.

We consider an indoor VLP system as illustrated in Figure 1, which is deployed
in a room with the size of $L\times W\times H$, where $N_{\mathrm{t}}$ LEDs
are attached at the ceiling as positioning anchors, with the position
coordinates of
$\mathbf{a}_{i}=[x_{i}^{\mathrm{a}},y_{i}^{\mathrm{a}},z_{i}^{\mathrm{a}}]$.
Without loss of generality, the optical wireless channel of both the light-of-
sight (LOS) link and the non-light-of-sight link with once reflection on the
walls are considered (Komine and Nakagawa, 2004), and the optical channel gain
between the $i$-th LED and the $j$-th UE located at
$\mathbf{c}_{j}=[x_{j},y_{j},z_{j}]$ is given by

(1) $H_{i,j}=H_{i,j}^{\mathrm{LOS}}+H_{i,j}^{\mathrm{NLOS}},$

where $H_{i,j}^{\mathrm{LOS}}$ and $H_{i,j}^{\mathrm{NLOS}}$ represents the
channel direct current (DC) gains of the LOS and NLOS links, respectively.

Assuming that the LEDs follow a Lambertian emission pattern, the channel gain
of the LOS link can be expressed as

(2)
$H_{i,j}^{\mathrm{LOS}}=\left\\{\begin{array}[]{cc}\frac{(m+1)A_{\mathrm{PD}}\cos^{m}\left(\varphi_{i,j}\right)\cos\left(\theta_{i,j}\right)T_{\mathrm{s}}g\left(\theta_{i,j}\right)}{2\pi
d_{i,j}^{2}},&0\leq\theta_{i,j}\leq\theta_{\mathrm{FOV}},\\\
0,&\theta_{i,j}>\theta_{\mathrm{FOV}},\end{array}\right.$

where $d_{i,j}=\left\|\mathbf{c}_{j}-\mathbf{a}_{i}\right\|_{2}$ is the Euclid
distance between the $i$-th LED and the $j$-th UE. $m$ is the Lambertian
order, $A_{\mathrm{PD}}$ is the effective physical area of the PD.
$T_{\mathrm{s}}$ and $g\left(\theta_{i,j}\right)$ are the gain of the optical
filter and the optical concentrator of the PD, respectively.
$\theta_{\mathrm{FOV}}$ is the field-of-view (FOV) angle of the PD.
$\varphi_{i,j}$ and $\theta_{i,j}$ are the angles of irradiance and incidence
with respect to the normal direction, respectively. In addition, the
Lambertian order is given by $m=-\ln 2/\ln(\cos(\varphi_{1/2}))$, where
$\varphi_{1/2}$ is the half-power angle of the LED. And the gain of the
optical concentrator $g\left(\theta_{i,j}\right)$ can be expressed as

(3)
$g\left(\theta_{i,j}\right)=\left\\{\begin{array}[]{cc}\frac{n^{2}}{\sin^{2}(\theta_{\mathrm{FOV}})},&0\leq\theta_{i,j}\leq\theta_{\mathrm{FOV}},\\\
0,&\theta_{i,j}>\theta_{\mathrm{FOV}},\end{array}\right.$

where $n$ denote the refractive index of the optical concentrator.

The NLOS link is modeled as a diffusion channel, where the wall is segmented
into many small diffuse reflectors with the area of $A_{\mathrm{ref}}$ and the
average diffuse reflectance of $\bar{\rho}$. The NLOS channel gain can be
approximately calculated as

(4)
$H_{i,j}^{\mathrm{NLOS}}=\left\\{\begin{array}[]{cc}\sum_{\mathrm{Wall}}\bar{\rho}h_{i,j}^{\mathrm{in}}h_{i,j}^{\mathrm{ref}},&0\leq\theta_{i,j}\leq\theta_{\mathrm{FOV}},\\\
0,&\theta_{i,j}>\theta_{\mathrm{FOV}}.\end{array}\right.$

Specifically, $h_{i,j}^{\mathrm{in}}$ and $h_{i,j}^{\mathrm{ref}}$ are the
channel gains of incidence and reflection light given by

(5)
$h_{i,j}^{\mathrm{in}}=\frac{A_{\mathrm{ref}}(m+1)}{2\pi(d_{i,j}^{\mathrm{in}})^{2}}\cos^{m}\left(\varphi_{i,j}^{\mathrm{in}}\right)\cos\left(\theta_{i,j}^{\mathrm{in}}\right),$
(6)
$h_{i,j}^{\mathrm{ref}}=\frac{A_{\mathrm{PD}}}{2\pi(d_{i,j}^{\mathrm{ref}})^{2}}\cos\left(\varphi_{i,j}^{\mathrm{ref}}\right)\cos\left(\theta_{i,j}^{\mathrm{ref}}\right),$

where $d_{i,j}^{\mathrm{in}}$ represents the distance between the $i$-th LED
and the diffuse reflector, and $d_{i,j}^{\mathrm{ref}}$ represents the
distance between the diffuse reflector and the $j$-th PD.
$\varphi_{i,j}^{\mathrm{in}}$ and $\theta_{i,j}^{\mathrm{in}}$ are the angle
of irradiance and incidence of the incident light, respectively.
$\varphi_{i,j}^{\mathrm{ref}}$ and $\theta_{i,j}^{\mathrm{ref}}$ are the angle
of irradiance and incidence of the reflective light, respectively.

Thus, the received electronic power from the $j$-th LED is given by

(7)
$P_{i,j}=\left[R_{\mathrm{p}}P_{i}^{\mathrm{t}}\left(H_{i,j}^{\mathrm{LOS}}+H_{i,j}^{\mathrm{NLOS}}\right)\right]^{2}+\sigma_{\mathrm{n}}^{2},$

where $P_{i}^{\mathrm{t}}$ is the optical power emitted by the $i$-th LED,
$R_{\mathrm{p}}$ is the responsivity of the PD.
$\sigma_{\mathrm{n}}^{2}=\sigma_{\mathrm{shot}}^{2}+\sigma_{\mathrm{thermal}}^{2}$
is the noise variance composed of two components of the shot noise and the
thermal noise, which are given by

(8)
$\sigma_{\mathrm{shot}}^{2}=2qR_{\mathrm{p}}P_{i,j}B+2qI_{\mathrm{bg}}I_{2}B,$
(9) $\sigma_{\mathrm{thermal}}^{2}=\frac{8\pi kT_{\mathrm{K}}}{G}\eta
AI_{2}B^{2}+\frac{16\pi^{2}kT_{\mathrm{K}}\Gamma}{g_{\mathrm{m}}}\eta^{2}A^{2}I_{3}B^{3},$

where $q$ is the electronic charge, and $k$ is the Boltzmann constant.
$I_{\mathrm{bg}}$ is the background current of the ambient light. $B$ is the
noise bandwidth. $T_{\mathrm{K}}$ is the absolute temperature. $G$ is the
open-loop voltage gain. $\eta$ is the fixed capacitance of the PD. $\Gamma$
and $g_{\mathrm{m}}$ are the channel noise factor and the transconductance of
the field effect transistors, respectively. $I_{2}$ and $I_{3}$ denote the
noise bandwidth factors, which are defined as 0.562 and 0.0868, respectively
(Tang et al., 1996). Finally, the received power vector of the $j$-th UE can
be modeled as $\mathbf{p}_{j}=[P_{1,j},P_{2,j},\cdots,P_{N_{\mathrm{t}},j}]$.

###  2.2. Problem Statement

As mentioned above, the received power vector $\mathbf{p}_{j}$ is dependent on
the channel gain $H_{i,j}$, which is controlled by the position-related
parameters. Thus, the correlation between $\mathbf{c}_{j}$ and
$\mathbf{p}_{j}$ can be exploited to estimate the location coordinate of the
UE. Then, the positioning task can be mathematically formulated as

(10)
$\min_{\mathbf{w},\mathcal{M}(\cdot)}\mathbb{E}_{\mathbf{c}_{j}\in\mathbb{C}}\left\|\hat{\mathbf{c}}_{j}-\mathbf{c}_{j}\right\|_{2}^{2},$

where $\hat{\mathbf{c}}_{j}=\mathcal{M}(\mathbf{p}_{j},\mathbf{w})$ denotes
the estimated position coordinate obtained by a model $\mathcal{M}(\cdot)$
with the learnable parameters of $\mathbf{w}$. $\mathbb{C}$ is the set of the
positions in the whole indoor environment lighten by the LEDs. To solve this
problem, on the one hand, an efficient model $\mathcal{M}(\cdot)$ is needed to
extract the position-related features and map them into the label space of
location coordinates. On the other hand, it is also important to find a set of
optimal parameters $\mathbf{w}^{*}$ adaptive to the indoor environment, which
is more challenging in nonstationary environments, because of the highly time-
varying channel condition caused by the environmental changes, such as ambient
light interference, obstruction movement, and aging hardware. To improve the
positioning accuracy and robustness in nonstationary environments, in this
paper, an FL-based cooperative VLP scheme is proposed, which is introduced in
detail in Section 3.

##  3\. Proposed Cooperative VLP Scheme

###  3.1. Cooperative VLP Framework

Figure 2. The cooperative VLP framework based on FL.

To deal with the nonstationary environments, in this paper, a cooperative VLP
framework is designed, as illustrated in Figure 2. A typical indoor scenario
in the hospital is considered as an example, where $N_{\mathrm{t}}$ LEDs are
deployed as access points of both VLC and VLP, and $N_{\mathrm{r}}$ UEs
equipped with abundant computing resources, e.g., mobile phones or intelligent
medical devices, are carried by users and move in the building. Using the
uplink wireless transmission, such as WiFi and Bluetooth, the UEs can transmit
the data of local model weights to the server. Moreover, each UE can store a
local dataset and a local model for position estimation, which is feasible for
off-the-shelf mobile devices (Ignatov et al., 2018). In the proposed
cooperative VLP framework, UEs participate in the collection of training data
and construct their local datasets $\mathbf{\Omega}_{j}$. Exploiting the local
datasets, the UEs cooperatively train a global model applicable for the whole
indoor environment using the FL-based algorithm, which is elaborated in the
following subsection.

###  3.2. FL-Based VLP Algorithm

In this section, we introduce the proposed FL-based VLP algorithm, starting
from the multi-user data collection, and then describing the structure of the
VCPosNet and the FL-based training and inference procedures. The block diagram
of the proposed FL-based VLP algorithm is illustrated in Figure 3.

Figure 3. The block diagram of the proposed FL-based VLP algorithm.

####  3.2.1. Multi-User Sensing and Position Data Collection

To ensure the samples in local datasets can represent the real-time indoor
environment, UEs collect the training samples and update the local datasets
over time. Specifically, in a sensing interval, the $j$-th UE measures the
received signal power vector
$\mathbf{p}_{j}^{d}=[P_{1,j}^{d},P_{2,j}^{d},\cdots,P_{N_{\mathrm{t}},j}^{d}]$
transmitted from the $N_{\mathrm{t}}$ LEDs, and labels them with the
corresponding position coordinate
$\mathbf{c}_{j}^{d}=[x_{j}^{d},y_{j}^{d},z_{j}^{d}]$. As users change their
locations, the received signal power vector $\mathbf{p}_{j}^{d}$ corresponding
to different locations are collected by their UEs. Thus, each UE can construct
a local dataset
$\mathbf{\Omega}_{j}=\\{\mathbf{p}_{j}^{d},\mathbf{c}_{j}^{d}\\}_{d=1}^{D_{j}}$
consisting of $D_{j}$ samples. Unlike the most learning-based VLP methods
using a one-shot site survey, in the proposed scheme, data collection is
implemented over time. Then, the local dataset is accordingly updated with the
new samples, to make it suitable for the real-time environment.

####  3.2.2. VCPosNet Structure for Federated Learning

To better extract the location-related features contained in the received
power vector $\mathbf{p}_{j}$ and improve the positioning accuracy, we devise
a neural network architecture called VCPosNet in this paper.

In the proposed VCPosNet, firstly, a convolution layer (CL) with
$N_{\mathrm{f}}$ kernel filters is employed, of which the output is applied
with the rectified linear units (ReLU) activation function. Then, to prevent
possible overfitting because of accumulated parameters, some neurons are
dropped out with a probability in the training stage. After this, the
resultant data is normalized to ensure that the output data has zero mean and
unit standard deviation, to accelerate network training and prevent gradient
vanishing. Subsequently, the output data is fed into the long short-term
memory (LSTM) layer consisting of $N_{\mathrm{u}}$ LSTM units, whose output is
fed into a Tanh activation function, and also followed by a dropout and a
normalization layer. Finally, a fully-connected layer is connected with the
final output of the neural network to generate the estimated position
coordinate $\hat{\mathbf{c}}_{j}$.

####  3.2.3. FL-Based Continuous Model Training and Inference

To obtain a neural network model adaptive to the nonstationary environments,
in this paper, FL is applied to train the VCPosNet in a cooperative manner.
Specifically, the training process is composed of the local training processes
at the UEs and the federated aggregation of the global model at the server.
The detailed procedures are performed as follows and illustrated in Algorithm
1.

Algorithm 1 FL-Based VLP Algorithm (Training Stage)

0: 1) Maximum number of local epochs $E$2) Minibatch size $B$3) Local datasets
$\boldsymbol{\Omega}_{1},\boldsymbol{\Omega}_{2},\cdots,\boldsymbol{\Omega}_{N_{\mathrm{r}}}$

1: Initialize the global weights $\mathbf{w}_{0}$

2: for each communication round $t=1,2,\cdots$ do

3: At UEs: (parallelly performed by each UE $j$, $\forall j=1,\cdots
N_{\mathrm{r}}$)

4: Initialize the local weights $\mathbf{w}^{(j)}\leftarrow\mathbf{w}_{t}$

5: for each local epoch $e=0,\cdots,E-1$ do

6: Sample a mini-batch $\boldsymbol{\Omega}_{j}^{e}$ from the local dataset
$\boldsymbol{\Omega}_{j}$

7: Compute the loss function
$\mathcal{L}(\mathbf{w}^{(j)},\boldsymbol{\Omega}_{j}^{e})$ by (11) and update
the local weights $\mathbf{w}^{(j)}$ using back-propagation

8: end for

9: Submit the local weights $\mathbf{w}_{t+1}^{(j)}=\mathbf{w}^{(j)}$ and
dataset size $D_{j}$ to server via the uplink transmission

10:

11: At the Server:

12: Aggregate the weights
$\mathbf{w}_{t+1}^{(1)},\cdots,\mathbf{w}_{t+1}^{(N_{\mathrm{r}})}$ by (12)
and broadcast $\mathbf{w}_{t+1}$ to all the $N_{\rm{r}}$ UEs via VLC downlink
transmission

13: end for

At each UE, the local dataset $\mathbf{\Omega}_{j}$ is used to train the
weights of VCPosNet to obtain a local model containing the environment
information of local areas. Mean square error (MSE) is utilized as the loss
function, which is given by

(11)
$\mathcal{L}\left(\boldsymbol{\Omega}_{j},\mathbf{w}\right)=\frac{1}{D_{j}}\sum_{d=1}^{D_{j}}\left\|\hat{\mathbf{c}}_{j}^{d}-\mathbf{c}_{j}^{d}\right\|_{2}^{2},$

where $\hat{\mathbf{c}}_{j}^{d}=\mathcal{M}(\mathbf{p}_{j}^{d},\mathbf{w})$
represents the output of the VCPosNet with the input of $\mathbf{p}_{j}^{d}$
and the model weights of $\mathbf{w}$. In each local epoch, the local weights
$\mathbf{w}^{(j)}$ are updated using stochastic gradient descent and back-
propagation to minimize the loss function using mini-batches
$\boldsymbol{\Omega}_{j}^{e}$ of size $B$ sampled from
$\boldsymbol{\Omega}_{j}$. After $E$ local epochs, the $j$-th UE uploads the
learned local weights $\mathbf{w}_{t+1}^{(j)}$ and the dataset size $D_{j}$ to
the server.

At the server, by applying a weighted average proportional, the local weights
$\mathbf{w}_{t+1}^{(1)},\mathbf{w}_{t+1}^{(2)},\cdots,\mathbf{w}_{t+1}^{(N_{\mathrm{r}})}$
received from UEs are aggregated into the global weight $\mathbf{w}_{t+1}$,
which is calculated by

(12)
$\mathbf{w}_{t+1}=\sum_{j=1}^{N_{\mathrm{r}}}\frac{D_{j}}{D}\mathbf{w}_{t+1}^{(j)},$

where $D=D_{1}+D_{2}+\cdots+D_{N_{\mathrm{r}}}$ is the total number of
training samples owned by all the UEs. Then, the global weights
$\mathbf{w}_{t+1}$ are broadcasted to all the UEs via the VLC downlink
transmission.

Afterwards, the $j$-th UE updates its local weight $\mathbf{w}^{(j)}$ with the
received global weights $\mathbf{w}_{t+1}$, and then continues the local
training processes in the next communication round. With the increase of
communication rounds, a temporally generalized position estimation model
covering the whole LED-illuminated indoor environment is well trained with the
cooperation of multiple users. Moreover, once the indoor environment changes,
the global model weights will converge to the optimal values again using
updated local datasets in a few communication rounds. Thus, the global model
will adapt to the current channel condition rapidly, making it robust to the
environmental changes.

Since network training is continuously conducted, we consider that the
inference process may happen in any communication round. Once the $j$-th UE
requests for the instantaneous position coordinate, the $j$-th UE can measure
the received power vector $\mathbf{p}_{j}$ with PD, then input it into the
local model with the latest global weights $\mathbf{w}_{t}$ to estimate the
real-time location coordinate $\hat{\mathbf{c}}_{j}$.

With the aid of the FL framework and the cooperation of multiple users, an
adaptive neural network with robust generalization and adaptation ability can
be obtained, which is verified by the simulation results in Section 4. In
addition, as the training datasets containing massive private data are kept
locally, the privacy of users is effectively protected through physical
isolation.

##  4\. Simulation Result and Discussion

Figure 4. Positioning error with respect to communication rounds.

Figure 5. CDF with respect to the positioning error.

Figure 6. Positioning error in nonstationary environments.

In this section, extensive simulations have been conducted to evaluate the
performance of the proposed scheme. We consider an indoor scenario in a room
with the size of $L\times W\times H=5\times 5\times 3\mathrm{m}^{2}$, where
$N_{\mathrm{t}}=16$ LEDs are vertically mounted at the inner vertexes of
$5\times 5$ latticed ceiling. The transmitted optical power of the LED is set
as $P_{i}^{\mathrm{t}}=1\mathrm{W}$, and the half-power angle is
$\varphi_{1/2}=60^{\circ}$. The Lambertian order is $m=1$. In addition, the
average reflectance of the wall is set as $\bar{\rho}=0.7$. As to the PD, the
FOV angle is set as $\theta_{\mathrm{FOV}}=90^{\circ}$. The effective area is
set as $A_{\mathrm{PD}}=1\mathrm{cm}^{2}$, and the responsivity is set as
$R_{\mathrm{p}}=0.6\mathrm{A/W}$. Moreover, the optical filter gain is set as
$T_{\mathrm{s}}=1$, and the optical concentrator refractive index is set as
$n=1$. The background current is set as $I_{\mathrm{bg}}=740\mu\mathrm{A}$,
and the noise bandwidth is set as $B=5\mathrm{MHz}$. Other parameters are set
following literature (Komine and Nakagawa, 2004). We assume that
$N_{\mathrm{r}}=10$ UEs participate in cooperative model training. Each UE
owns a local dataset with the size of $D_{\mathrm{r}}=900$, and updates the
local dataset with newly collected data every 10 communication rounds. In the
training process, the number of local epochs is set as $E=5$, and the
minibatch size is set as $B=128$.

The positioning error of the proposed FL-based VLP scheme with respect to
communication rounds is shown in Figure 6, compared with benchmark schemes of
MLP, CNN, recurrent neural network (RNN), and LSTM trained using the FL
framework. It is shown by the results that, the positioning error of the
proposed scheme decreases rapidly over time with the cooperation of multiple
UEs, and reaches centimeter-level accuracy after only two communication
rounds, which verifies the superior convergence capability of the proposed FL-
based VLP scheme. After 200 communication rounds, the positioning error of the
proposed scheme reaches 3.38 cm, which outperforms the benchmark schemes by
6.52 cm, 3.61 cm, 4.67 cm, and 4.37 cm, respectively.

Subsequently, the cumulative distribution function (CDF) with respect to the
positioning error of the proposed scheme is reported in Figure 6, and compared
with benchmark methods, including the deep-learning-based methods of MLP, CNN,
RNN, and LSTM, and machine-learning-based methods of KNN and support vector
machine (SVM). The deep-learning-based methods are trained with the FL
framework. The machine-learning-based methods are usually implemented relying
on the whole dataset and thus are trained in a centralized manner. It is noted
that the CDF of the proposed scheme grows rapidly, and then reaches the value
of 0.96 at the positioning error of 10 cm, which means that the positioning
accuracy can be controlled within the centimeter-level with a probability of
96%. Moreover, it can be observed that, the proposed scheme can achieve a
higher CDF at the different positioning errors than benchmarks, which verifies
the proposed VCPosNet can extract the position-related knowledge effectively
and further improve the positioning accuracy.

To further investigate the generalization capability of the proposed method in
nonstationary environments, we consider three dynamically changing scenarios,
i.e., some parameters of the indoor environment change over time. The details
of the three scenarios are as follows. (A) Time-varying ambient light
intensity. The background current of the ambient light is initially set as
$I_{\mathrm{bg}}=0$ before the first communication round, and is increased by
$50\mu\mathrm{A}$ after every communication round. (B) Hardware failure or
obstruction blocking. Some LED signals are unavailable to the UEs, i.e., the
received signal power from an LED suddenly changes to zero before the first
communication round. (C) Device aging. The LED emitting power degrades
gradually over time, which is modeled as
$P_{i}^{\mathrm{t}}=P^{\mathrm{init}}\exp(-\beta t)$, where the initial power
$P^{\mathrm{init}}$ is set as 1 W, and the decay coefficient $\beta$ is set as
$-\ln(0.8)/100$.

The positioning errors of the proposed FL-based VLP scheme in nonstationary
environments are illustrated in Figure 6, compared with the centralized-
learning-based scheme based on the one-shot site survey. For a fair
comparison, the initial model weights of both the proposed scheme and the
centralized-learning-based scheme are set as the pre-trained model weights
before the first communication round. It is shown by the results that, the
positioning accuracy of the centralized-learning-based scheme degrades rapidly
with the aggravation of environmental changes. However, in the proposed FL-
based VLP framework, the global model weights can adapt to the time-varying
channel conditions with the cooperation of multiple users, and the positioning
error can be kept at a satisfactory level, which verifies the generalization
capability in nonstationary environments.

##  5\. Conclusion

In this paper, a privacy-preserving cooperative VLP scheme based on FL is
proposed, where multiple users cooperatively train a global model adaptive to
nonstationary environments. Additionally, The CVPosNet network structure is
presented to extract the position-related features and improve the positioning
performance. Simulation results show that the proposed scheme outperforms the
benchmark schemes in positioning accuracy and generalization performance,
especially in the nonstationary environment.

## References

  * (1)
  * Arfaoui et al. (2021) Mohamed Amine Arfaoui, Mohammad Dehghani Soltani, Iman Tavakkolnia, Ali Ghrayeb, Chadi M. Assi, Majid Safari, and Harald Haas. 2021\.  Invoking Deep Learning for Joint Estimation of Indoor LiFi User Position and Orientation.  _IEEE Journal on Selected Areas in Communications_ 39, 9 (2021), 2890–2905.
  * Ignatov et al. (2018) Andrey Ignatov, Radu Timofte, William Chou, Ke Wang, Max Wu, Tim Hartley, and Luc Van Gool. 2018.  AI Benchmark: Running Deep Neural Networks on Android Smartphones. In _Proceedings of the European Conference on Computer Vision (ECCV) Workshops_. 288–314.
  * Keskin et al. (2018) Musa Furkan Keskin, Ahmet Dundar Sezer, and Sinan Gezici. 2018\.  Localization via Visible Light Systems.  _Proceedings of the IEEE_ 106, 6 (2018), 1063–1088.
  * Komine and Nakagawa (2004) Toshihiko Komine and Masao Nakagawa. 2004.  Fundamental analysis for visible-light communication system using LED lights.  _IEEE Transactions on Consumer Electronics_ 50, 1 (2004), 100–107.
  * Konečný et al. (2016) Jakub Konečný, H. Brendan McMahan, Daniel Ramage, and Peter Richtárik. 2016.  Federated Optimization: Distributed Machine Learning for On-Device Intelligence.  _CoRR_ (2016), 1–38.  arXiv:1610.02527
  * Lin and Zhang (2020) Xintong Lin and Lin Zhang. 2020.  Intelligent and Practical Deep Learning Aided Positioning Design for Visible Light Communication Receivers.  _IEEE Communications Letters_ 24, 3 (2020), 577–580.
  * Liu et al. (2015) Sicong Liu, Fang Yang, Chao Zhang, and Jian Song. 2015\.  Narrowband Interference Mitigation Based on Compressive Sensing for OFDM Systems.  _IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences_ E98.A, 3 (2015), 870–873.
  * Luo et al. (2017) Junhai Luo, Liying Fan, and Husheng Li. 2017.  Indoor Positioning Systems Based on Visible Light Communication: State of the Art.  _IEEE Communications Surveys & Tutorials_ 19, 4 (2017), 2871–2893.
  * Matheus et al. (2019) Luiz Eduardo Mendes Matheus, Alex Borges Vieira, Luiz F. M. Vieira, Marcos A. M. Vieira, and Omprakash Gnawali. 2019.  Visible Light Communication: Concepts, Applications and Challenges.  _IEEE Communications Surveys & Tutorials_ 21, 4 (2019), 3204–3237.
  * Tang et al. (1996) Andrew P. Tang, Joseph M. Kahn, and Keang-Po Ho. 1996.  Wireless infrared communication links using multi-beam transmitters and imaging receivers. In _Proceedings of International Conference on Communications_. 180–186.
  * Van et al. (2017) Manh The Van, Nguyen Van Tuan, Tran The Son, Hoa Le-Minh, and Andrew Burton. 2017.  Weighted k-nearest neighbour model for indoor VLC positioning.  _IET Communications_ 11, 6 (2017), 864–871.
  * Yang et al. (2016) Fang Yang, Junnan Gao, and Sicong Liu. 2016.  Priori Aided Compressed Sensing-Based Clipping Noise Cancellation for ACO-OFDM Systems.  _IEEE Photonics Technology Letters_ 28, 19 (2016), 2082–2085.
  * Yang et al. (2017) Fang Yang, Junnan Gao, Sicong Liu, and Jian Song. 2017\.  Clipping Noise Elimination for OFDM Systems by Compressed Sensing With Partially Aware Support.  _IEEE Transactions on Broadcasting_ 63, 1 (2017), 103–110.
  * Yang et al. (2019a) Fang Yang, Junnan Gao, and Yuanpeng Liu. 2019a.  Indoor visible light positioning system based on cooperative localization.  _Optical Engineering_ 58, 1 (2019), 1–6.
  * Yang et al. (2019b) Qiang Yang, Yang Liu, Tianjian Chen, and Yongxin Tong. 2019b.  Federated Machine Learning: Concept and Applications.  _ACM Transactions on Intelligent Systems and Technology_ 10, 2 (2019), 1–19.
  * Zafari et al. (2019) Faheem Zafari, Athanasios Gkelias, and Kin K. Leung. 2019\.  A Survey of Indoor Localization Systems and Technologies.  _IEEE Communications Surveys & Tutorials_ 21, 3 (2019), 2568–2599.
  * Zhuang et al. (2018) Yuan Zhuang, Luchi Hua, Longning Qi, Jun Yang, Pan Cao, Yue Cao, Yongpeng Wu, John Thompson, and Harald Haas. 2018.  A Survey of Positioning Systems Using Visible LED Lights.  _IEEE Communications Surveys & Tutorials_ 20, 3 (2018), 1963–1988.