Limits on astrophysical antineutrinos with the KamLAND experiment

We report on a search for electron antineutrinos ($\bar{\nu}_e$) from astrophysical sources in the neutrino energy range 8.3 to 30.8 MeV with the KamLAND detector. In an exposure of 6.72 kton-year of the liquid scintillator, we observe 18 candidate events via the inverse beta decay reaction. Although there is a large background uncertainty from neutral current atmospheric neutrino interactions, we find no significant excess over background model predictions. Assuming several supernova relic neutrino spectra, we give upper flux limits of 60--110 cm$^{-2}$ s$^{-1}$ (90% CL) in the analysis range and present a model-independent flux. We also set limits on the annihilation rates for light dark matter pairs to neutrino pairs. These data improves on the upper probability limit of $^{8}$B solar neutrinos converting into $\bar{\nu}_e$'s, $P_{\nu_e \rightarrow \bar{\nu}_e}<3.5\times10^{-5}$ (90% CL) assuming an undistorted $\bar{\nu}_e$ shape. This corresponds to a solar $\bar{\nu}_e$ flux of 60 cm$^{-2}$ s$^{-1}$ (90% CL) in the analysis energy range.


INTRODUCTION
Underground liquid-scintillator neutrino detectors observe geo neutrinos, solar neutrinos, and reactor neutrinos below 10 MeV energy, in addition to the atmospheric neutrinos peak at O(GeV) range. However, other astrophysical neutrino sources also exist in our universe: from supernova explosions to hypothetical dark-matter annihilation neutrinos. The valley of the neutrino energy spectrum between the end of reactor neutrinos and the onset of atmospheric neutrinos can be used to search for astrophysical neutrinos. We present a search for astrophysical neutrinos in the neutrino energy between 8.3 and 30.8 MeV, focusing on electron antineutrinos (ν e ) from the Sun, past supernovae, and dark matter annihilation.
The Sun is the dominant source of astrophysical neutrinos, and various neutrino detectors have observed solar electron neutrinos (ν e ) (Abe et al. 2011;Aharmim et al. 2013;Abe et al. 2016;Agostini et al. 2020). As discussed in Malaney et al. (1990), antineutrinos are also produced in the Sun in comparatively small amounts from beta decays of 40 K, 232 Th, and 238 U, and they has not been observed yet. Solar neutrinos can be converted to antineutrinos with combined processes of the Mikheyev-Smirnov-Wolfenstein (MSW) effect (Smirnov 2005) and Resonant Spin Flavor Precession (RSFP) (Lim & Marciano 1988;Akhmedov 1988) as discussed in Akhmedov & Pulido (2003) and Díaz et al. (2009). This happens in two-step processes: The RSFP is a neutrino helicity resonance transition similar to the MSW effect in the Sun via the neutrino magnetic moment (µ). A simple RSFP model is excluded due to the large neutrino magnetic moment required, µ > 10 −10 µ B , already excluded by experiments (Beda et al. 2013;Agostini et al. 2017). The µ B is the Bohr magneton. The combined RSFP+MSW model is still allowed. The conversion probability is expressed as P (ν e →ν e ) 1.8 × 10 −10 sin 2 2θ 12 × µ 10 −12 µ B B T (0.05R ) 10 kG where θ 12 is the neutrino mixing angle in the Pontecorvo-Maki-Nakagawa-Sakata matrix, B T is the transverse solar magnetic field in the region of neutrino production, R is the solar radius, µ is the neutrino magnetic moment. Experimentally, the conversion probability for solar 8 B neutrinos was studied by KamLAND (Gando et al. 2012), Borexino (Agostini et al. 2021), and Super-K (Abe et al. 2020).
A supernova explosion is one of the largest neutrino burst events in our universe. Supernova neutrinos from SN1987A, which occurred on 1987 February 23rd in the Large Magellanic Cloud, were detected by the water cherenkov detectors: KamiokaNDE (Hirata et al. 1987(Hirata et al. , 1988 and IMB (Bionta et al. 1987;Bratton et al. 1988), and the Baksan scintillation detector (Alexeyev et al. 1988). A future nearby supernova explosion could reveal detailed information on the explosion mechanism. At the same time, neutrinos from all the past supernovae are still traveling in our universe. These are called supernova relic neutrinos (SRN), and they provide the diffuse supernova neutrino flux. The SRN energy spectrum and associated detection rates have been discussed in various models (Kaplinghat et al. 2000;Horiuchi et al. 2009;Nakazato et al. 2013Nakazato et al. , 2015. The most stringent experimentalν e flux upper limit is given by Super-K (Abe et al. 2021a), but no significant signal observation has been made yet. KamLAND is able to perform a comparable search forν e at around 10 MeV. The higher neutron tagging efficiency should give an advantage over Super-K searching at the lower energy region.
Neutrinos can also be produced in the annihilation of dark matter particles. In case of the existence of an MeV-scale light dark matter particle, its self-annihilation process might produce neutrino pairs (χχ → νν) at MeV energies. Assuming a model of MeV-scale dark matter annihilation in our galactic halo (Palomares-Ruiz & Pascoli 2008), theν e flux from dark-matter self-annihilation is given by where m χ is the dark matter mass, σ A v is the averaged self-annihilation cross section times the relative velocity of the annihilating particles, J ave is the angular-averaged intensity over the whole Milky Way, R sc = 8.5 kpc is the distance of the Sun to the galactic center, and ρ 0 = 0.3 GeV cm −3 is the local dark matter density. Here, a factor 1/3 is assumed for the branching ratio to the three flavors of neutrinos. This process is also discussed in Klop & Ando (2018) and Argüelles et al. (2019).
In this work, we show the results for two benchmark cases of J ave = 1.3 and 5.0 (Palomares-Ruiz & Pascoli 2008). In this paper, after describing the KamLAND experiment (Section 2), we describe the search forν e with an energy range of 8.3-30.8 MeV (Section 3). The main backgrounds are discussed in Section 4; these are reactor neutrinos, accidental backgrounds, spallation products, fast neutrons, and atmospheric neutrinos. The data is interpreted in Section 5 and we present the conversion probability of solar 8 B antineutrinos, the SRN flux, and the dark-matter self-annihilation cross section. We summarize those results in Section 6.

KAMLAND DETECTOR
The KamLAND experiment uses ultra pure liquid scintillator to detectν e via the inverse beta-decay (IBD) reaction. The detector is located 1000 m underground, underneath Mt. Ikenoyama in Kamioka Japan, corresponding to 2700 m water equivalent. The cosmic muon flux is suppressed by ∼ O(10 5 ) relative to sea level. Figure 1 shows a schematic view of the detector. KamLAND consists of a 18-m diameter stainless-steel sphere tank (Inner Detector, ID) and a cylindrical vessel of water-cherenkov muon veto (outer detector, OD) surrounding the ID. A 13-m diameter EVOH/nylon balloon (outerballoon) holds 1-kton of liquid scintillator at the center of the ID. Photon sensors, 1325 17-inch and 554 20-inch photomultiplier tubes (PMTs), sensitive to scintillation light are bolted to the inside of the stainless-steel tank. Details of the detector design are described in Suzuki (2014).
KamLAND data acquisition started in March 2002, and this study uses all data sets up to July 2020. The KamLAND-Zen 400 phase of the project operated with a 154 cm-radius nylon balloon (inner-balloon) at the center of KamLAND filled with xenon-loaded liquid scintillator from August 2011 to December 2015 (Gando et al. 2016). The detector was further upgraded in May 2018 when KamLAND-Zen 800 started with a new 1.9 m-radius inner-balloon (Gando 2020;Gando et al. 2021). In this study, the inner-balloon volume is vetoed to avoid possible background contamination. The OD system was refurbished in 2016, when the 225 PMTs were replaced by 140 new PMTs including 47 higher quantum efficiency PMTs (Ozaki & Shirai 2017).
The interaction vertex and energy deposition are reconstructed using the measured PMT charge and timing information. At low energies, the detector was calibrated using various radioactive sources: 60 Co, 68 Ge, 203 Hg, 65 Zn, 241 Am 9 Be, 137 Cs, and 210 Po 13 C. Above 10 MeV, the energy response is calibrated using spallation-products of 12 B (τ = 29.1 ms, Q = 13.4 MeV) and 12 N (τ = 15.9 ms, Q = 17.3 MeV).
The position-dependent energy calibration and fiducial volume determination uncertainty were studied with calibration sources positioned throughout the outer balloon volume (Berger et al. 2009). The reconstructed energy and interaction vertex resolution were determined to be 6.4%/ E (MeV) and ∼ 12 cm/ E (MeV) (Gando et al. 2013), respectively. Daily stability measurements were performed using 2.2 MeV gamma rays emitted from spallation-induced neutron capture on protons (Abe et al. 2010) and spallation 12 B events. The total estimated uncertainty including the time variation of the energy scale, linearity, and uniformity is within ±2.0% for this data set.
The primary radioactive backgrounds in the liquid scintillator are (5.0±0.2)×10 −18 g g −1 of 238 U and (1.8 ± 0.1) × 10 −17 g g −1 of 232 Th (Gando et al. 2015). These radioactive contaminants are negligibly small relative to other backgrounds in this study, such as muon spallation products and atmospheric neutrinos, and are therefore ignored.

ELECTRON ANTINEUTRINO SELECTION
Electron antineutrinos are detected in KamLAND via the IBD reaction (ν e + p → e + + n) with a 1.8 MeV neutrino energy threshold. IBD candidate events are selected by the delayed coincidence (DC) method: scintillation light from the positron and its annihilation gamma-rays is the prompt event, followed by a 2.2 MeV (4.9 MeV) gamma-ray from neutron capture on a proton (carbon-12) after a mean capture time of 207.5 ± 2.8 µs (Abe et al. 2010) (delayed event). The incident neutrino energy (E ν ) is computed from the reconstructed prompt energy (E prompt ), E ν E prompt + 0.8 MeV + E n , where E n is the average neutron kinetic energy of O(10 keV).
The DC selection criteria between the prompt and delayed events use the prompt energy, delayed energy (E delayed ), spatial distribution (∆R), and time difference (∆T ); they are 7.5 < E prompt < 30 MeV, 1.8 < E delayed < 2.6 MeV or 4.4 < E delayed < 5.6 MeV, ∆R < 160 cm, and 0.5 < ∆T < 1000 µs, respectively. The two delayed-energy selection criteria correspond to a 2.2-MeV capture gamma-ray on a proton and a 4.9-MeV capture gamma-ray on a carbon-12. The timing difference between the prompt and delayed events are required from the 207.5-µs of neutron capture time. The spatial correlation selection is optimized from diffusion length of a thermalized neutron and delayed gamma ray. The selection efficiency is evaluated with Monte Carlo (MC) simulations described in next paragraph. In this energy range, one of the primary backgrounds is a fast neutron (described in Sec. 4), mostly in close proximity to the ID vessel. In order to suppress this contamination, we select a 550 cm-radius spherical fiducial volume from the center of KamLAND, corresponding to a total number of target protons of N p = (4.6 ± 0.1) × 10 31 .
During the KamLAND-Zen 400/800 phases, the inner-balloon regions are vetoed for the delayed event in order to avoid background contamination from the xenon-loaded liquid scintillator, inner balloon body, and suspending ropes. The inner-balloon regions are cut from the analysis: a 250cm-radius spherical volume centered in the detector and a 250-cm-radius vertical cylindrical volume in the upper-half of the detector. For the above DC selection, the IBD detection efficiency IBD is estimated through MC simulations with uniformly generated neutrino events and is determined to be IBD = 92%(73%) for without(with) inner-balloon cut.
The data period from May 2002 to July 2020 totals 4528.5 days of livetime. We find 21 DC pairs after DC selection, 3 of them have multiple delayed events following the prompt event and are excluded from the final sample as they are likely due to the fast neutron backgrounds and/or atmospheric neutrino interactions. Figure 2 and Figure 3 present the electron antineutrino candidate distributions. Possible backgrounds in this analysis come from reactor antineutrinos, the accidental coincidence of events, spallation products, fast neutron, and atmospheric neutrinos. Radioactive backgrounds and reactor neutrinos were studied during the reactor-and geo-neutrino measurements at Kam-LAND (Gando et al. 2011a(Gando et al. , 2013, while previousν e analyses in the O(10 − 10 3 ) MeV energy range (Gando et al. 2012;Asakura et al. 2015a,b;Abe et al. 2021b) showed that fast neutron and atmospheric neutrinos are the most challenging backgrounds above ∼ 10 MeV.

Reactor antineutrinos
The location of the KamLAND detector is surrounded by 56 Japanese nuclear power reactors. The reactor neutrino flux comes primarily from the beta decay of neutron-rich fragments produced in the fission of 4 isotopes: 235 U, 238 U, 239 Pu, and 241 Pu. For each reactor, the appropriate operational records including thermal power generation, fuel burn-up, shutdowns, and fuel reload schedule were used to calculate the fission rates. With the reactor operation data, we have measured reactor antineutrinos between a few to several MeV (Gando et al. 2011b(Gando et al. , 2013. However, there are no measured reactor neutrino spectra in the present analysis energy range. Hence, we use the reactor neutrino spectra assuming polynomial functions based on the results from the ILL experiment results (Huber 2011;Mueller et al. 2011). We find that the background of this analysis contribution from reactor neutrinos becomes negligibly small above 10 MeV. In the range, it has a large spectrum shape uncertainty of ∼ 50% from the Huber/Mueller spectrum model. The number of reactor neutrino backgrounds is estimated to be 1.4 ± 0.6 including KamLAND-detector related uncertainties. The expected number of events from the extrapolated reactor spectrum is consistent with the Daya Bay results ) at neutrino energies of [8.125,12 MeV], the highest energy bin in the Daya Bay analysis.

Accidental coincidence
Two uncorrelated events may accidentally pass through the DC selection. Predominantly, uncorrelated long-lived spallation isotopes or radioactive decays could produce a mimic prompt event, and radioactive decays such as 214 Bi and 208 Tl beta/gamma-rays possibly become a 2.2 MeV of mimic delayed event. To estimate the random coincidence background, events were selected with appropriate prompt and delayed energies but in an off-time window of 0.2-1.2 s after the prompt event. This offtime window is 10 3 times larger than the antineutrino selection, providing a high statistics estimate. The expected number of accidental coincidence background events is (7.3 ± 1.0) × 10 −2 .

Spallation products
Cosmic muons induce various spallation products in KamLAND (Abe et al. 2010). Short-lived spallation products are rejected by a 2 ms whole volume veto. Some longer-lived products are a potential background in this study. As a primary spallation cut, we apply a 2 ms whole volume veto for all muons in KamLAND, which has a muon rate of ∼ 0.34 Hz. Muons in the ID are identified when more than O(10 4 ) photons are detected by the PMTs. We identify muon events in the OD when the number of OD hits exceeds 5 or 9 hits, before and after the OD refurbishment campaign, respectively (Ozaki & Shirai 2017).
The previous analysis of KamLAND data (Gando et al. 2012) showed that the 9 Li (τ = 257.2 ms, Q = 13.6 MeV) spallation product is a challenging background. In order to reduce this background contamination, we improved the spallation veto introducing a likelihood-ratio based muon shower tagging in addition to the primary 2 ms veto. Using a similar idea employed in Super-K analysis (Bays et al. 2012), we evaluate a probability density function for spallation-like events taking into consideration the spatial and timing correlation of the muon track and charge deposition. Due to low statistics of 9 Li in KamLAND data with a production rate of 2.8 ± 0.2 kton −1 day −1 (Abe et al. 2010), it is difficult to directly estimate the correlation between the muon track and 9 Li. Hence, instead of the 9 Li, we use 12 B data, whose production rate is ∼ 20 times larger. Figure 4 (a) shows the closest track distance distribution (dL) between the muon track and the spallation production point for 12 B, 9 Li, and neutrons, based on a Monte-Carlo (MC) study with FLUKA (version 2011.2x.8.patch) (Böhlen et al. 2014;Ferrari et al. 2005) and propagated through KamLAND with Geant4 (version 4.9.6 patch-04) (Agostinelli et al. 2003;Allison et al. 2006;Allison et al. 2016). The spread of the 9 Li distance distribution is narrower than 12 B. This means that a tight spallation cut on the 12 B data can be used to put an upper limit on the remaining 9 Li background. The difference of muon charge deposition among the spallation products was small enough. Figure 4 (b) shows the data-driven correlation between muon charge deposition per track length (dE/dX) and the distance of spallation products from the muon track (dL). Muons depositing a large charge in the detector induce spallation products even far from the muon track. Considering the lifetime of 9 Li, a 2 s veto is sufficient to reject the background but the detector livetime becomes too small. Here we define a likelihood-ratio parameter depending on the dE/dX, dL, and the time difference from the muon event to the subsequent event in order to optimize the rejection of spallation events. To optimize the likelihood-ratio threshold while maximizing detector livetime and minimizing the spallation-cut inefficiency, we define the figure of merit (FOM) as follows (Punzi 2003): where livetime is the detector livetime ratio, N spall. (non spall.) are the expected number of spallation (non-spallation) events without a spallation veto, cut is the spallation cut efficiency depending on the likelihood-ratio threshold, δ 2 cut is the cut efficiency uncertainty, and δ 2 stat. is the statistical uncertainty on the expected number of events. Maximizing the FOM with the condition that the 9 Li spallation cut inefficiency becomes zero consistent within the range determined from 12 B, we optimize the likelihoodratio threshold. This muon-shower based likelihood-ratio spallation cut allows livetime = 79% of detector livetime on average and gives a spallation cut inefficiency of (0.2 ± 0.5)%, that means 99.8% of spallation background reduction. This gives a spallation background of 1.4 ± 3.6 in this analysis energy range. Figure 4. (a) The closest distance distribution from the muon track to the spallation products production points of 12 B, 9 Li, and neutron from FLUKA and Geant4 simulations. Only muons leaving more than 10 6 photo-electrons were selected in this analysis. Dashed lines correspond to the cumulative ratio (righthanded scale). 99% of these spallation products are within 300 cm from the track. (b) Data-driven correlation between muon charge deposition (dE/dX) and the closest distance of spallation product from the muon track (dL). The color bar shows the likelihood-ratio value.

Fast neutrons
The fast neutrons background comes from outside of the detector and is induced by cosmic muons in the surrounding rock and water. Neutron scattering on protons or carbon nuclei in the liquid scintillator can mimic a prompt event. After that, the neutrons are thermalized and captured on a proton or carbon, creating the delayed event. The 2 ms veto for OD tagged muons rejects the majority of this background but some events remain due to OD inefficiency.
This background was evaluated with a Geant4-based MC, which included a detailed description of the KamLAND geometry (KLG4sim). Neutron interactions were treated with the QGSP BIC HP physics list, while muon-nucleus interactions were activated using G4EmExtraPhysics. The cosmic muon directional distributions were implemented from the KamLAND spallation simulation study (Abe et al. 2010) which used a topological map of Mt. Ikenoyama (Geographical SurveyInstitute of Japan 1997) and the MUSIC simulation tool (Antonioli et al. 1997). The simulated detector response in KLG4sim was tuned with various calibration data.
An equivalent of 8313 livetime days were simulated in KLG4sim. In the case of a muon going through the ID producing a lot of scintillation emission, the muon and associated neutrons were vetoed by the 2 ms whole volume veto. Figure 5 (a) shows the reconstructed fast neutron position distribution for OD-tagged MC events. For comparison, OD-tagged fast neutron events in the data are also shown in Figure 5 (a). The fast neutron selection used the IBD selection described in Section 3 except for the OD tagging and ∆T > 10 µs selection avoiding decay-electron contribution. While the data has a slightly broader spread due to the difficulty of vertex reconstruction around the boundary between liquid scintillator and buffer oil, the radial distribution in the fiducial volume are consistent between data and simulation. The fast neutron radial distribution f (R) is assumed to be f (R) ∝ exp (R/λ) as a function of the distance R from the detector center, where λ = (50.9 ± 3.0) cm. We use this non-uniform position distribution in the fit to data to evaluate the fast neutron background (see Section 5). The fast neutron energy spectra in MC and data are shown in Figure 5   The remaining fast neutrons generate events with a few OD hits below the OD detection threshold which are therefore not tagged in the OD (OD-untagged). From the KLG4sim simulation, the estimated number of fast neutron background within the fiducial volume is 6.8 ± 6.8 DC pairs by scaling the detector livetime. The uncertainty is conservatively estimated to be 100%, considering the poorly known production yield of neutrons in the rock. Events producing multiple neutrons via the 12 C(n, 2n) 11 C reaction are expected to give 3.4 ± 3.4 DC pairs.

Atmospheric neutrinos
Atmospheric neutrinos produce the dominant background in this analysis via charged current (CC) and neutral current (NC) interactions. In the CC interactions, various neutron emission modes are possible backgrounds. Since the neutrino-nucleus interaction cross section for carbon is at least one order of magnitude smaller than for protons (Kim & Cheoun 2009), the background comes mainly from proton interactions. IBD from atmospheric electron antineutrinos is not a dominant contribution in the CC background because the mean energy is higher than our analysis energy range and small flux in this range. Based on the atmospheric neutrino flux (Honda et al. 2007) and measured crosssection by MiniBooNE (Athar et al. 2007), the dominant process isν µ + p → µ + + n. In the KamLAND detector, the muon decay is observed as muon scintillation, muon decay, and neutron capture. The two-prompts (muon scintillation + its decay) and one-delayed (neutron capture) DC events are vetoed with ∼ 78% efficiency (Gando et al. 2012). Its inefficiency contributes to the background. The number of CC backgrounds is estimated to be 1.1 ± 0.3.
To estimate the NC background, we took into account the atmospheric neutrino flux (Honda et al. 2007), cross sections (Ahrens et al. 1987), the neutron binding energies in carbon for the P-shell (18.7 MeV) and the S-shell (41.7 MeV) configurations and the corresponding shell populations, and de-excitation models (Kamyshkov & Kolbe 2003). The NC interaction is given by ν(ν) + 12 C → ν(ν) + n + 11 C + γ. Most of the outgoing neutrons have a kinetic energy of less than 200 MeV and they scatter on protons resulting in a visible energy of typically less than 100 MeV. The details of the background signatures and estimations are described in Gando et al. (2012). The resulting expected number of NC interactions in this data set is 20.6 ± 5.9, where the uncertainty comes from the atmospheric neutrino flux and the cross section which are combined to provide ∼ 30% in total.
For comparison, we also estimate the atmospheric neutrino background with NEUT (version 5.4.0.1) (Hayato 2009;Hayato & Pickering 2021). The interaction models are summarized in Table 1. We use the Honda et al. (2015) model for the atmospheric neutrino flux including the matter oscillation effect implemented in the Prob3++ (Wendell 2012). The de-excitation model for oxygen is incorporated in NEUT, but that for carbon is missing. After the final state interaction in NEUT, the outgoing particles were introduced into KLG4sim. The response of this simulation was compared to KamLAND data in the 200 MeV-1.5 GeV energy range, outside of the fast neutron background. Although there are large uncertainties from the fast neutrons in the 30-200 MeV energy range, data and MC are consistent within the errors. Below 100 MeV, the NC quasi-elastic scattering (NCQE) is dominant. The NC two-particle-two-hole (NC 2p2h) interaction will also contribute, but assuming that the ratio of the NC 2p2h to the NCQE cross sections is similar to the corresponding the CC ratio, roughly 5-10% (Nieves et al. 2011), the contribution is estimated to be small compared to the NCQE contribution and its large uncertainties.
For the DC energy selection of 7.5-30 MeV in NEUT based estimation, the remaining CC and NC backgrounds are estimated to be 0.9 +0.3 −0.4 and 16.5 +5.1 −4.5 , respectively. The NEUT background estimate are smaller, but are consistent within the uncertainties. The neutron multiplicity in the NC reaction may play a role in the models and affect the background estimate due to the DC requirement of selecting a single neutron capture. The energy spectrum shape is concrete on the de-excitation models of carbon and the proton-scattering by neutrons in the numerical calculation, but this simulationbased estimation does not include the de-excitation. Therefore we took the numerically calculated spectrum (Gando et al. 2012) and estimated the number of NC backgrounds to be 20.6 ± 5.9 in this analysis. We treat the number of NC background events as a free parameter in this analysis, independent the number of backgrounds from the estimation models, and use the energy spectrum to constrain.

ANALYSIS AND RESULTS
Our search for astrophysical electron antineutrino signals fitted the energy spectra and radial distributions in data to the estimated backgrounds. The fast neutron background contributes with a large uncertainty but is mostly concentrated at the outer radius, while the other backgrounds and neutrino candidates have a uniform distribution in the detector. The atmospheric neutrino NC interaction is the primary background in this analysis. We used the following χ 2 to fit the number of atmospheric NC backgrounds and the number of astrophysical neutrinos: with, where N observed is the number of observed IBD candidates, N astro.ν is the number of astrophysical neutrino events, N NC is the number of atmospheric neutrino NC background events, N BG i (i = 1, 2, ...5) represent the number of the other background contributions (see Table 2). The statistical uncertainty σ stat. is the square root of the total number of expected events. In the shape χ 2 term (χ 2 shape ), R is the radius, E is the energy, f j (R) is the normalized radius distribution, and g j (E) the normalized energy spectrum for each contribution j where j = 1, 2, ...7 correspond to the astrophysical neutrino signal, atmospheric NC background, and the other 5 background contributions. We use them as an unbinned log-likelihood fit to test χ 2 shape . Only fast neutrons have an exponential radius distribution f i (R), the other contribution are uniform. We integrate the energy and radius over 7.5-30 MeV and 0-550 cm, respectively. The penalty term (χ 2 penalty ) is computed from the systematic uncertainties (δ k ): energy spectrum shape uncertainty, radial distribution uncertainty, detector efficiency uncertainty, and energy scale uncertainty. In the background term (χ 2 BG ), N BG i is the number of the i-th background events, N expected BG i is the expected number of the i-th background component, and δ 2 BG i is its associated uncertainty.

Solar electron antineutrino
Assuming an unoscillated 8 B neutrino flux of 5.94 × 10 6 cm −2 s −1 (Pena-Garay & Serenelli 2008), the region allowed by the fit is shown in Figure 6 and summarized in Table 2. The best fit values for the ν e →ν e conversion probability and NC events are 0.0 and 7.5 ± 3.4, respectively. The number of atmospheric neutrino NC interactions is smaller than the estimate but model 2σ and data 2σ bands overlap. This value is also consistent with the NEUT simulation result within 1σ. Figure 7 shows the energy and radial distributions for best-fit backgrounds and the upper limit for solar 8 Bν e with 90% confidence level (CL). All residual values are within ±2σ region. The obtained upper limit on the conversion probability is 3.5 × 10 −5 at 90% CL, corresponding to a 60 cm −2 s −1 solar 8 Bν e flux limit above 8.3 MeV of neutrino energy (containing 30% of the solar 8 B neutrino flux). In a comparable case of using a measurement 8 B neutrino flux of 5.25 × 10 6 cm −2 s −1 (Aharmim et al. 2013), the upper limit on the conversion probability becomes 3.9 × 10 −5 at 90% CL. This result improves on the previous KamLAND study (Gando et al. 2012) and is the most stringent upper limit to date.
From the upper limit on the conversion probability and Equation (3), we also obtain the upper limit on the neutrino magnetic moment (µ) and the transverse solar magnetic field (B T ) in the region of neutrino production: µ < 4.9 × 10 −10 µ B 10 kG using 34 • for the mixing angle θ 12 (Gando et al. 2011b). This bound is weaker than the most stringent upper limit of 0.28 × 10 −10 µ B from the solar neutrino spectrum measurement by Borexino (Agostini et al. 2017). Figure 6. The results and allowed regions for the solar ν e →ν e conversion probability and the number of atmospheric neutrino NC interactions. Color contours correspond to 1σ (red), 90% (orange), 2σ (green), and 3σ (blue). The best-fit conversion probability and NC events are 0 and 7.5, respectively (black dot). The horizontal hatched region represent the expected number of NC events with 1σ uncertainty. Top and right panels are 1-dimensional ∆χ 2 distributions for conversion probability (CP) and number of atmospheric neutrino NC interaction, respectively. The upper limit on the conversion probability is 3.5×10 −5 at 90% CL.

Supernova relic neutrinos
To search for SRNs, we fit with different theoretical models that produceν e emission: the Kaplinghat model (Kaplinghat et al. 2000) (Kaplinghat+00), the Horiuchi model in the case of 6 MeV effective temperature (Horiuchi et al. 2009) (Horiuchi+09, 6 MeV), the Nakazato maximum model in the case of inverted mass ordering (Nakazato+15 max, IH), and the Nakazato minimum model in the case of normal mass ordering (Nakazato+15 min, NH) (Nakazato et al. 2013(Nakazato et al. , 2015. From the χ 2 defined in Equation (6), we find no significant excess of SRNs with any of the models. As an example of the fitting result with the Nakazato+15 (max, IH) model, the best-fit value for the number of SRNs is 0 events while the number of NC backgrounds is 7.5 events. This result is consistent with the calculated number of SRN events of 0.4 in KamLAND. The 90% CL upper limit on the number of  events is 9.3. The upper flux limit is calculated to be 108 cm −2 s −1 from where F 90 and N 90 are the upper limits on flux and the number of events, respectively. The dF dE M and dN dE M are the theoretical differential flux and spectrum, respectively, for the SRN models. This 90% CL flux upper limit is still much higher than the expected flux of 5.1 cm −2 s −1 . Table 3 shows a summary of the fit results for each theoretical model and corresponding upper limit. For all tested models, the best fit number of SRN is 0, and NC background is 7.5. However, the reported upper limit changes for each model due to differences in the underlying theoretical energy spectrum.  (Nakazato et al. 2013(Nakazato et al. , 2015 9.3 108 5.1 Nakazato+15 (min, NH) (Nakazato et al. 2013(Nakazato et al. , 2015 8.9 105 2.2 Note-F 90 and N 90 are the 90% CL upper limits of flux and number of events, respectively. The expected flux is integrated over our analysis energy range E prompt = [7.5, 30 MeV].

Model independent flux
We also present model-independent upper limits on theν e flux assuming monochromatic neutrino energies. The flux upper limits (φ 90 ) are calculated with where N 90 is the 90% CL upper limit on the number ofν e in a 1 MeV wide bin using the Feldman & Cousins (1998) approach, N p is the number of target protons, σ is the IBD reaction cross section, IBD is the detector efficiency, and T is the detector livetime. Figure 8 shows the resulting electron antineutrino flux in comparison with results from Borexino (Agostini et al. 2021), Super-K (Bays et al. 2012;Zhang et al. 2015;Abe et al. 2021a), and various theoretical SRN models (Kaplinghat et al. 2000;Horiuchi et al. 2009;Nakazato et al. 2013Nakazato et al. , 2015. While our results do not yet exclude SRN models, they provide the strictest flux limits for E ν = [8.3, 13.3 MeV]. Table 4 shows a summary of the flux upper limits per a bin.

Dark matter self annihilation
Theν e flux upper limit per energy bin can be translated to a dark matter self-annihilation cross section limit (Palomares-Ruiz & Pascoli 2008). From Equation (4), we obtain an upper limit of σ A v < (1-11) × 10 −24 cm −3 s −1 (90% CL) for the benchmark case of J ave = 1.3 (see Figure 9). This result is the most stringent constraint on the self-annihilation cross section for m χ < 15 MeV.

SUMMARY
We searched for astrophysicalν e in the neutrino energy range 8.3 to 30.8 MeV with 4528.5 livetime days of KamLAND data. No significant excess was found over the expected backgrounds. We presented the strictest upper limit on the conversion probability of solar 8 B neutrinos to antineutrinos, 3.5 × 10 −5 (at 90% CL). Assuming various model predictions, the upper limit on the SRN flux translates to 60-110 cm −2 s −1 . We also give the strictest upper limit on the model independent flux below 13.3 MeV but this limit is still an order of magnitude larger than SRN model predictions. The upper limits on the dark matter self-annihilation cross-section to neutrino pairs are the most stringent for dark matter particle masses below 15 MeV. Our results for the model-independent flux limit (Table 4) can set limits on various astrophysicalν e 's, for instance, neutrinos from sterile neutrino decay (Hostert & Pospelov 2020) and primordial black hole dynamics (Dasgupta et al. 2020;Wang et al. 2021;Calabrese et al. 2021).
Further background suppression is necessary to improve the solar 8 Bν e and SRN sensitivity. A future neutrino detector at a deep underground site will suppress the spallation background (Anderson et al. 2019;Guo et al. 2021). A larger distance to nuclear power plants will reduce the reactor neutrino component. More detailed measurements of the high-energy reactor neutrino spectrum are necessary, including the end point. The background arising from atmospheric neutrinos is the most challenging. Pulse shape discrimination (PSD) may reduce this contribution (Li et al. 2016) in a future large neutrino detector such as JUNO . Although the fast scintillation decay time of the current KamLAND liquid scintillator cocktail and significant re-emission, PSD can be improved by the detector upgrades and requires excellent timing resolution for PMTs. The KamLAND2 detector upgrade program intends to use a linear-alkyl-benzene based liquid scintillator which would realize the PSD due to a slower scintillation decay time compared to the current KamLAND liquid scintillator (Asakura et al. 2015c;Obara et al. 2019;Kamei 2020;Nakamura et al. 2020;Takeuchi & Kawada 2020).
The KamLAND experiment is supported by JSPS KAKENHI Grants 18J10498, 19H05102, and 19H05803; the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan; Netherlands Organization for Scientific Research (NWO); and under the U.S. Department of Energy (DOE) Contract No. DE-AC02-05CH11231, the National Science Foundation (NSF) No. NSF-1806440, NSF-2012964, as well as other DOE and NSF grants to individual institutions. The Kamioka Mining and Smelting Company has provided service for activities in the mine. We acknowledge the support of NII for SINET4. We also thank Y. Hayato for advising our atmospheric neutrino simulation with NEUT. This work is partly supported by the Graduate Program on Physics for the Universe (GP-PU).