Search for Solar Flare Neutrinos with the KamLAND Detector

Solar flares are the largest explosions in the solar system, releasing energy between 10²⁸–10³³ erg in only tens of minutes (Schrijver et al. 2012). The mechanism of solar flares can be described as the rapid conversion of magnetic energy to thermal and kinetic energy of charged particles by reconnection of the magnetic field on the solar surface (Parker 1957). Observations of electromagnetic signals, ranging from radio waves to γ-rays at 100 MeV, and neutrons emitted during solar flares contribute to the current understanding of this phenomenon (Benz 2008).

In the standard flare model, solar flares accelerate protons to more than 300 MeV and then nuclear reactions of accelerated protons generate pions in the solar atmosphere (Hudson & Ryan 1995). The decay of these pions produces high-energy (>70 MeV) γ-rays and MeV-scale neutrinos. Thus, neutrino production is expected in the standard solar flare model, and the properties of these solar flare neutrinos depend on the initial accelerated proton spectrum and flux (Kocharov et al. 1991).

In recent decades, neutrino emission models from solar flares have been developed, and such models inform the feasibility of detecting solar flare neutrinos. Fargion (2004) predicted that detection of neutrinos from a large solar flare (>10³² erg) was feasible with Super-Kamiokande and IceCube. Recent updates, however, predict no possibility of detecting solar flare neutrinos even with Hyper-Kamiokande (Takeishi et al. 2013). Another study (de Wasseige 2016) predicts 398–770 cm⁻² neutrino fluence at Earth in the 10–100 MeV range, which corresponds to ≪1 electron scattering in KamLAND. From these recent studies (Takeishi et al. 2013; de Wasseige 2016), it is clear that MeV neutrino observation from a single flare is hardly feasible. However, by searching for a statistical excess in coincidence with a large number of solar flares, it may be possible to detect solar flare neutrinos. Such a detection can provide an additional probe to understand the particle acceleration on the solar surface.

There have been several efforts to experimentally search for solar flare neutrinos. The Homestake experiment reported a small excess of events correlated with a large solar flare in 1991 (Davis 1994). On the other hand, KAMIOKANDE II and LSD observed no excess of events associated with different solar flares (Hirata et al. 1990; Aglietta et al. 1991). The Sudbury Neutrino Observatory (SNO) has performed a coincidence search with 842 solar flares measured from radiation from 3 keV to 17 MeV with the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) and found no correlations (Aharmim et al. 2014). In 2019, an analysis by Borexino improved the upper limits on neutrino fluence and excluded the Homestake parameter space (Agostini et al. 2021). In this analysis, the Borexino collaboration assumed that the neutrino flux is proportional to X-ray intensity and used 472 M- and X-class solar flares selected from the Geostationary Operational Environmental Satellite (GOES) database. The aforementioned studies are sensitive to neutrinos in the 1–100 MeV range. Recently, IceCube reported the first search for GeV-scale neutrinos related to intense γ-ray solar flares and constrained some of the parameter space associated with theoretical predictions for the neutrino flux (Abbasi et al. 2021).

In this paper, we present a search for solar flare neutrinos using the KamLAND data taken from 2002 March to 2019 September, which includes solar cycles 23 and 24. KamLAND is a 1 kton liquid-scintillator detector that is sensitive to neutrinos in the energy range between 1 MeV and a few GeV. However, in this study, we focus on 1–35 MeV neutrinos. For experimental studies of solar flare neutrinos, flare selection and the time window for coincidence studies are important. We discuss these in Section 2. Section 3 provides an overview of the KamLAND detector and the two detection channels for our solar flare neutrino search. The scheme of the coincidence analysis is presented in Sections 4 and 5. The analysis results are converted to fluence upper limits in Section 6.

In a solar flare, neutrinos are generated from charged pion decay. Neutral pion decay emits 70–100 MeV γ-rays. It is natural to identify solar flares and set the coincidence time window from the γ-ray measurements (de Wasseige 2016). IceCube applied this strategy to the γ-ray data taken by the Fermi-LAT satellite (Abbasi et al. 2021).

However, as the Fermi-LAT satellite was launched in 2008, solar γ-ray burst observations are not available for the 23rd solar cycle, including the largest (class X28) flare on record that occurred on 2003 November 4. For this reason, we apply another strategy to identify solar flares and the timing of particle acceleration in solar flares.

Hard X-ray is an alternative channel to identify flares and set the time windows. Hard X-ray emission is generated from the bremsstrahlung of nonthermal electrons accelerated to relativistic velocities by a solar flare. Shih et al. (2009) reported that there is a close proportionality between the line γ-ray and hard X-ray fluence from solar flares. The existence of line γ-rays is indirect evidence of hadronic interactions in a solar flare, which is one of potential sources of neutrino emission.

The light curve of solar flare X-rays differentiated by time is expected to be similar to the light curve of hard X-rays through the Neupert effect (Neupert 1968; Dennis & Zarro 1993). The Neupert effect is an experimentally known effect that the derivative of soft X-ray light curves from solar flares tends to have the same timing response as microwave emission. It can also be applied in the case of hard X-rays instead of microwave emission. Thus, we can find the time window for the solar flare neutrino search using the differential soft X-ray light curves from GOES satellites.

The advantages of the GOES soft X-ray profile compared to the RHESSI hard X-ray/γ-ray profile or the Fermi-LAT γ-ray profile are the length of the observation and the availability of stable data. With the method mentioned above, we can use the most abundant data set of the solar flare since 1975, which covers the 23rd and 24th solar cycles, and set the flare time window even if the hard X-ray and γ-ray observations are not available during the solar flare. This method is suggested and validated in Okamoto et al. (2020).

Based on Okamoto et al. (2020), we determine the flare time window as follows: (i) calculate the differential X-ray light curve, (ii) search for the peak of the differential curves, and (iii) define the time window starting from the nearest zero coefficient before the peak and ending at the nearest zero coefficient after the peak. The red curve in Figure 1 is one of the examples of our flare time window. The duration time of this example is 1143 s.

Zoom In Zoom Out Reset image size

We obtain the flare list from the GOES X-ray database at the National Oceanic and Atmospheric Administration. After the X- and M-class selection, which was also used in the Borexino analysis (Agostini et al. 2021), there were 1342 flares with a total X-ray intensity of 639.3 × 10⁻⁴ W m⁻² from 2002 March to 2019 September. For the coincidence analysis with the KamLAND data, all time windows were required to be in a period of operation in which the live time to running time ratio of the detector was more than 95%. The KamLAND live time is defined as the integrated period of time that the detector was sensitive to neutrinos and includes corrections for calibration periods, detector maintenance, daily run switch, etc. Applying these requirements, we found 614 solar flares remained. The distributions of the duration and intensity of these flares are shown in Figures 2 and 3, respectively. The average length of the 614 time windows is 1028 s. The duration time described in Figure 1 is almost the mean value. The integrated intensity is 303.0 × 10⁻⁴ W m⁻², which is 25 times larger than the flare coincident with the Homestake excess and 1.7 times larger than the flares used in the Borexino analysis (Agostini et al. 2021).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The KamLAND detector is a large-volume neutrino detector, which is located approximately 1 km underground under Mt. Ikenoyama in Kamioka, Japan. KamLAND consists of an outer water-Cerenkov detector and an inner scintillation detector. The water-filled outer detector (OD), housed in a 10 m radius × 20 m high cylindrical vessel, provides shielding from external γ-ray backgrounds and an active muon counter. The OD was instrumented with 225 20-inch Photo Multiplier Tubes (PMTs) before a refurbishment in 2016 and 140 20-inch PMTs after the refurbishment (Ozaki & Shirai 2017). The inner detector is a 9 m radius stainless-steel spherical tank with 1325 17-inch PMTs and 554 20-inch PMTs mounted on the inner surface. The main volume of the inner detector is a 1 kton liquid scintillator supported by a 6.5 m radius nylon/EVOH balloon installed at the center of the stainless-steel tank. This nylon/EVOH balloon is called the outer balloon. Outside the outer balloon is filled with nonscintillating buffer oil. Another smaller nylon balloon for KamLAND-Zen is called the inner balloon and is described later. The details of the KamLAND detector are described in Suzuki (2014).

KamLAND began data taking in 2002 March. From 2011 August, KamLAND started the KamLAND-Zen phase to search for the neutrinoless double-beta decay of ¹³⁶Xe using a nylon balloon (inner balloon) installed at the center of the detector; this inner balloon is filled with a xenon-loaded liquid scintillator (Gando et al. 2016). During the initial phase, known as KamLAND-Zen 400, which ran from 2011 August to 2015 September, the inner-balloon radius was 1.5 m and the mass of xenon was about 400 kg. In 2018 May, the KamLAND-Zen experiment was upgraded to the so-called KamLAND-Zen 800 phase, with an enlarged inner balloon of radius 1.9 m and double the amount of xenon (about 800 kg) for a higher-sensitivity search (Gando 2020; Gando et al. 2021). For the KamLAND-Zen periods, the regions with a xenon-loaded scintillator were excluded from the effective volume for the neutrino search to suppress backgrounds from the xenon nuclei, nylon balloon, and supporting structures.

KamLAND has multiple reaction channels to detect neutrinos. We use the following two channels: neutrino–electron elastic scattering (ES), ν + e⁻ → ν + e⁻, and inverse-beta decay (IBD), ${\bar{\nu }}_{e}+p\to {e}^{+}+n$. ES is sensitive to all flavors of neutrinos, though the cross section depends on the neutrino flavor. This channel does not provide a measurement of the neutrino energy, though the energy of the scattered electron provides a lower bound. IBD is sensitive only to electron antineutrinos above 1.8 MeV. The IBD cross section is roughly 10 times larger than the ES cross section. In addition, the IBD signal has advantages to suppressing backgrounds thanks to a delayed coincidence measurement. The positron annihilates with an electron, emitting two 511 keV γ-rays. The positron and two γ-rays are observed as one event called the prompt event. The incident electron antineutrino energy, E_ν , can be reconstructed from the prompt scintillation as E_ν ≃ E_p + 0.8 MeV, where E_p is the energy of the prompt signal. With the mean capture time of about 207 μs, the neutron is captured on a proton (carbon) emitting a 2.2 (4.9)MeV γ-ray, which is called the delayed event. Exploiting the time-spatial correlation between the prompt and delayed events, we can observe electron-type antineutrinos in an almost background-free condition.

4.1. Basic Treatment of KamLAND Data

4.2. Selection Criteria for ES

Although there are some theoretical predictions of the spectrum of solar flare neutrinos (Kocharov et al. 1991; Fargion 2004), we conservatively assume a monochromatic spectrum for the solar flare neutrinos, like the gamma-ray burst neutrino analysis (Fukuda et al. 2002). For each assumed energy, E_ν , a lower-energy threshold (E_th) and analysis volume (V(r_fid)) were optimized to maximize the figure of merit (FoM), defined below. In this analysis, we used a spherical analysis volume, thus we optimized the analysis distance, r_fid, for the volume, V(r_fid). The FoM is defined as

$$\begin{eqnarray}&&\mathrm{FoM}=\displaystyle \frac{V({r}_{\mathrm{fid}})\times P({E}_{\mathrm{th}})}{\sqrt{B({r}_{\mathrm{fid}},{E}_{\mathrm{th}})}},\end{eqnarray} \tag{ 1 }$$

where P(E_th) is the probability that the energy of the ES electron exceeds E_th; B(r_fid, E_th) is the total number of background events in the flare-off time of that period with r < r_fid and ${E}_{\mathrm{th}}\lt {E}_{\mathrm{vis}}\lt {T}_{\max }$, where E_vis is the observed energy in the KamLAND detector; and T_max is the maximum kinetic energy of the recoil electron. This optimization was performed period by period. The detection efficiency, η^ES = V(r_fid)/V(600 cm) × P(E_th), resulting from the FoM optimization considering the detector energy-scale model is shown in Figure 4 as a function of the incident neutrino energy. The shape of η^ES(E_ν ) depends on the vertex distribution of external γ-ray backgrounds penetrating the tanks from the rock surrounding the detector.

Zoom In Zoom Out Reset image size

4.3. Background Estimation and χ² Studies

Around E_ν = 3 MeV, the background behavior changes. Below 3 MeV, there is a large number of backgrounds from radioactive decays in the balloons, PMTs, and the inner detector tank, such as ²⁰⁸Tl and ⁴⁰K. On the other hand, contributions from radioactivity are negligible above 3 MeV. Thus, we estimated the backgrounds above and below 3 MeV separately.

First, we describe the background determination below 3 MeV and flare coincidence analysis. Because the radioactive background rate was not sufficiently stable to estimate the rate in the solar flare time (on time) due to the liquid-scintillator convection on a timescale of hours, we estimated the accidental background in the following way. For the ith flare, 30 off-time windows were opened within a week before the flare. The duration of each off-time window was the same as of the on-time window. In each off-time window, the number of events with r < r_fid and ${E}_{\mathrm{th}}\lt {E}_{\mathrm{vis}}\lt {T}_{\max }$ was counted; these are shown as blue dots in Figure 5. We used the average (${N}_{i}^{\mathrm{off}}$) and standard deviation (σ_i ) of these off-time samples to estimate the expected number of background events with uncertainty for the associated on-time window. In Figure 5, ${N}_{i}^{\mathrm{off}}$ and σ_i are shown as a horizontal dashed line and a gray shaded region, respectively. Figure 5 is one example from the M1.8 flare in 2003 with E_ν = 1.0 MeV, where E_th = 0.4 MeV and r_fid = 600 cm. In this case, ${N}_{i}^{\mathrm{off}}$ and σ_i are 3567.2 and 80.0, respectively.

Zoom In Zoom Out Reset image size

The expected number of events in the on-time window with the solar flare signal of the ith flare is defined as ${n}_{i}\equiv {N}_{i}^{\mathrm{BG}}+{w}_{i}{\eta }^{\mathrm{ES}}{\alpha }^{\mathrm{ES}}{I}_{i};$ the first term represents the number of background events in the on-time window, and the second term corresponds to the number of signal events. The ${N}_{i}^{\mathrm{BG}}$ are assumed to follow a Gaussian distribution with a mean of ${N}_{i}^{\mathrm{off}}$ and a standard deviation of σ_i . In the second term, I_i is the flare intensity in units of [10⁻⁴ W m⁻²]; α^ES is a scale factor that connects the flare intensity and the number of ES in the 600 cm spherical volume, i.e., α^ES means how many electron scatterings occur in the 600 cm spherical volume by an X1 flare; η^ES is the detection efficiency described above; w_i is the detector live-time ratio in the ith on-time window. The observed number of events with r < r_fid and ${E}_{\mathrm{th}}\lt {E}_{\mathrm{vis}}\lt {T}_{\max }$ in the on-time window for the ith flare (${N}_{i}^{\mathrm{on}}$) should follow a Poisson distribution with a mean of n_i . In the case of Figure 5, ${N}_{i}^{\mathrm{on}}$ is 3525 and is shown as a red dot.

The χ² for all flares can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & 2\displaystyle \sum _{i\in \mathrm{flares}}[{N}_{i}^{\mathrm{on}}-{n}_{i}+{n}_{i}\mathrm{ln}({n}_{i}/{N}_{i}^{\mathrm{on}})]\\ & & +\displaystyle \sum _{i\in \mathrm{flares}}{\left(\displaystyle \frac{{N}_{i}^{\mathrm{BG}}-{N}_{i}^{\mathrm{off}}}{{\sigma }_{i}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 2 }$$

The second term in Equation (2) is a χ² penalty to account for the uncertainty on the background rate below 3 MeV, using σ_i as a conservative error. In Equation (2), α^ES and ${N}_{i}^{\mathrm{BG}}$ are free parameters, i.e., this χ² was minimized with respect to α^ES and ${N}_{i}^{\mathrm{BG}}(i=1,2,\cdots ,\,613)$.

Above 3 MeV, the background rate is small and stable. The mean number of events in the on-time window is ${n}_{i}=\langle {N}_{i}^{\mathrm{BG}}\rangle +{w}_{i}{\eta }^{\mathrm{ES}}{\alpha }^{\mathrm{ES}}{I}_{i}$, where $\langle {N}_{i}^{\mathrm{BG}}\rangle$ is the background rate averaged over the period scaled by the coincidence window duration. The χ² is modified to

$$\begin{eqnarray}&&{\chi }^{2}=2\sum _{i\in \mathrm{flares}}[{N}_{i}^{\mathrm{on}}-{n}_{i}+{n}_{i}\mathrm{ln}({n}_{i}/{N}_{i}^{\mathrm{on}})].\end{eqnarray} \tag{ 3 }$$

This χ² was minimized with respect to α.

From the χ² scan in our analysis range of 0.4–35 MeV for E_ν , the best-fit α^ES, ${\alpha }_{\mathrm{best}}^{\mathrm{ES}}$, was 0 for all assumed neutrino energies. The 90% confidence level (C.L.) upper limit on α^ES, ${\alpha }_{90}^{\mathrm{ES}}$, was estimated from ${\chi }^{2}({\alpha }_{\mathrm{best}}^{\mathrm{ES}})+2.7={\chi }^{2}({\alpha }_{90}^{\mathrm{ES}})$. The 90% confidence interval of α^ES is shown as a function of the assumed neutrino energy in Figure 6.

Zoom In Zoom Out Reset image size

We have a potential problem with our time window not matching the γ-ray emission time (Okamoto et al. 2020). To complement this, we also perform the coincidence analysis with a fixed time window. Another flare time window begins from the peak timing of the soft X-ray differential curve and runs for 1800 s. The complementary analysis shows the best-fit ${\alpha }_{\mathrm{best}}^{\mathrm{ES}}$ is consistent with zero within statistical errors for all assumed neutrino energies.

5.1. Selection Criteria for IBD

5.2. Background Estimation and χ² Studies

The IBD event rate is low and stable within each period because of the strong background reduction with the time–spatial correlation. Thus, a different χ² is defined as

$$\begin{eqnarray}&&{\chi }^{2}=2\sum _{p\in \mathrm{period}}[{N}_{p}^{\mathrm{on}}-{n}_{p}+{n}_{p}\mathrm{ln}({n}_{p}/{N}_{p}^{\mathrm{on}})],\end{eqnarray} \tag{ 4 }$$

where ${n}_{p}=\langle {N}_{p}^{\mathrm{off}}\rangle +{w}_{p}{\eta }^{\mathrm{IBD}}{\alpha }^{\mathrm{IBD}}{I}_{p}$ is the expected number of events in the cumulative flare time window, i.e., summed over all coincidence windows for the flares in our sample in the pth period. The ${N}_{p}^{\mathrm{off}}$ is the expected no-flare contribution, which is estimated from the IBD event rate in the pth period excluding the flare time window and scaled to the duration of the flare time window. ${N}_{p}^{\mathrm{on}}$ is the number of IBD events observed in the cumulative flare time window in the pth period. I_p is the cumulative X-ray intensity in the pth period. The parameter α^IBD is a scale factor that connects the flare intensity and the number of IBDs in the 600 cm spherical volume. In the IBD analysis, we used r_fid = 600 cm as the analysis distance and no energy binning to count events for the χ² study. The η^IBD indicates the detection efficiency for the electron antineutrinos via IBD and is computed with Monte Carlo simulation as shown in Figure 7. In the region below 4 MeV, the efficiencies are reduced due to larger accidental backgrounds that affect the likelihood selection. Because of the inner-balloon volume cuts during the KamLAND-Zen 400/800 phases as described in Section 3, the efficiencies in some periods are lower than those in other periods. Above about 4 MeV, the efficiencies converge to ∼77% for the inner-balloon cut periods and ∼94% for other periods. w_p is the detector live-time ratio.

Zoom In Zoom Out Reset image size

Assuming a monochromatic spectrum for the solar flare neutrinos, we varied E_ν from 1.8 MeV to 35 MeV and found the α^IBD that minimizes χ² for each assumed neutrino energies. The best-fit values of α^IBD, ${\alpha }_{\mathrm{best}}^{\mathrm{IBD}}$, and the 90% C.L. upper limits on α^IBD, ${\alpha }_{90}^{\mathrm{IBD}}$, were estimated with the same method in the ES analysis. The ${\alpha }_{\mathrm{best}}^{\mathrm{IBD}}$ was 0 for all assumed neutrino energies. The 90% confidence interval of α^IBD is shown in Figure 8.

Zoom In Zoom Out Reset image size

As in the ES analysis, we also performed a fixed-time window analysis and obtained ${\alpha }_{\mathrm{best}}^{\mathrm{IBD}}=0$ for all assumed neutrino energies.

Although there are some theoretical predictions of the spectrum of solar flare neutrinos (Kocharov et al. 1991; Fargion 2004), it has not been experimentally measured. Keeping the assumption of the monochromatic signal, we converted ${\alpha }_{90}^{\mathrm{ES}}$ into an upper limit on neutrino fluence, Φ^ES(E_ν ), as

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{\mathrm{ES}}({E}_{\nu })=\displaystyle \frac{{\alpha }_{90}^{\mathrm{ES}}({E}_{\nu })}{{N}_{{\rm{e}}}{\int }_{0}^{{T}_{\max }}\sigma ({E}_{\nu },{E}_{{\rm{e}}}){{dE}}_{{\rm{e}}}}\end{eqnarray} \tag{ 5 }$$

for the ES studies, where N_e is the number of electrons in the 6 m-radius spherical volume: 2.4 × 10³², E_e is the kinetic energy of the recoil electron, σ(E_ν , E_e) is the cross section of electron scattering with the incident neutrino of energy E_ν , and ${T}_{\max }$ is the maximum E_e. For the IBD studies, the upper limit on neutrino fluence, Φ^IBD(E_ν ), was obtained from

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{\mathrm{IBD}}({E}_{\nu })=\displaystyle \frac{{\alpha }_{90}^{\mathrm{IBD}}({E}_{\nu })}{{N}_{p}\sigma ({E}_{\nu })},\end{eqnarray} \tag{ 6 }$$

where N_p is the number of protons in the 6 m-radius spherical volume: (5.98 ± 0.13) × 10³¹, and σ(E_ν ) is the total cross section of the IBD from Strumia & Vissani (2003).

The fluence upper limit per flare is shown in Figure 9 with the assumption that all flares have the same neutrino luminosity, which is discussed in the SNO analysis (Aharmim et al. 2014). The allowed fluence region from the Homestake excess (Aharmim et al. 2014) is shown in the purple band. The upper limit by SNO (Aharmim et al. 2014) and Borexino (Agostini et al. 2021) are shown as the purple and green curves, respectively.

Zoom In Zoom Out Reset image size

Adopting another normalization that the solar flare neutrino luminosity is proportional to the X-ray intensity, the fluence upper limits were scaled to the Homestake flare's intensity (X12) as

$$\begin{eqnarray}&&{\rm{\Phi }}{({E}_{\nu })}^{\mathrm{scaled},\mathrm{ES}/\mathrm{IBD}}={{\rm{\Phi }}}^{\mathrm{ES}/\mathrm{IBD}}({E}_{\nu })\displaystyle \frac{12\times {10}^{-4}\,{\rm{W}}\,{{\rm{m}}}^{-2}}{303.3\times {10}^{-4}\,{\rm{W}}\,{{\rm{m}}}^{-2}}\end{eqnarray} \tag{ 7 }$$

and shown in Figure 10. The last term in Equation (7) represents the scaling factor from the flares analyzed in this work to the Homestake flare. The upper limit by Borexino (Agostini et al. 2021) is shown as a green curve after scaling to the X12 flare. Here, it is difficult to directly compare the results from KAMIOKANDE II (Hirata et al. 1990) and SNO (Aharmim et al. 2014) in this normalization because the flare catalog of those studies is different from this study, and we do not have sufficient X-ray measurements for the corresponding flares.

Zoom In Zoom Out Reset image size

In both normalizations, the 90% C.L. upper limits from this work exclude the entire region of the parameter space associated with the Homestake event excess for the large solar flare in 1991 and are the strictest upper limits on the solar flare neutrino fluence of all flavors in the neutrino energy range of 0.4–35 MeV.

We observe no evidence for neutrinos associated with solar flares in KamLAND. This work places the strictest upper limits on fluence normalized to the X12 flare with the assumption that the neutrino fluence is proportional to the X-ray intensity. At 20 MeV, the obtained 90% C.L limits are 8.4 × 10⁷ cm⁻² for electron antineutrinos and 3.0 × 10⁹ cm⁻² for electron neutrinos. The Homestake region is independently rejected by this result. To our knowledge, this is the first time the upper limit normalized to the flare intensity has been presented. We believe that this approach is useful to compare to results from other experiments and theoretical predictions.

The KamLAND experiment is supported by the JSPS KAKENHI grant 19H05803; the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan; Dutch Research Council (NWO); and under the US Department of Energy (DOE) contract No. DE-AC02-05CH11231, the National Science Foundation (NSF) Nos. NSF-1806440 and NSF-2012964, as well as other DOE and NSF grants to individual institutions. The Kamioka Mining and Smelting Company has provided services for activities in the mine. We acknowledge the support of N ii for SINET4. This work is partly supported by the Graduate Program on Physics for the Universe (GP-PU) and the Frontier Research Institute for Interdisciplinary Sciences, Tohoku University. A part of this study was carried out by using the computational resources of the Center for Integrated Data Science, Institute for Space-Earth Environmental Research, Nagoya University through the joint research program.