# Neutrino Oscillation Analysis of 217 Live Days of Daya Bay and 500 Live Days of RENO

Authors: Mario A. Acero; Alexis A. Aguilar-Arevalo; Dairo J. Polo-Toledo
Journal: Advances in High Energy Physics
 (2020)

---

## Abstract

We present a neutrino oscillation analysis of two particular data sets from the Daya Bay and RENO reactor neutrino experiments aiming to study the increase in precision in the oscillation parameterssin22θ13 and the effective mass splitting Δmee2 gained by combining two relatively simple to reproduce analyses available in the literature. For Daya Bay, the data from 217 days between December 2011 and July 2012 were used. For RENO, we used the data from 500 live days between August 2011 and January 2012. We reproduce reasonably well the results of the individual analyses, both rate-only and spectral, defining a suitable χ2 statistic for each case. Finally, we performed a combined spectral analysis and extract tighter constraints on the parameters, with an improved precision between 30 and 40% with respect to the individual analyses considered.


---

## Body

## 1. Introduction

Since their discovery in 1956 [1, 2], neutrinos have been under a heavy scrutiny by scientists trying to increase our knowledge about these abundant, exotic, and enigmatic particles. Neutrinos are neutral, spin−1/2, weakly interacting particles which are found to exist in three different flavors: electron neutrinos νe, muon neutrinos νμ, and tau neutrinos ντ. According to the SM, neutrinos are massless particles; however, a variety of experiments carried out over the past 50 years have shown that they undergo a quantum mechanical interference phenomenon, known as neutrino oscillation [3], through which the flavor of a neutrino changes while traveling from one point to another, implying that they must have nonzero masses. The discovery of neutrino oscillations was awarded the Nobel Prize in Physics in 2015.Within the standard theory of neutrino oscillations, a neutrino of a given flavor can be expressed as a superposition of three definite-mass neutrinosνkk=1,2,3 as
(1)να=∑kUαkνk,where Uαk are the elements of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix, which depend on the mixing angles θkj and a CP-violating phase δCP [4]. The PMNS matrix may also depend on two additional Majorana phases α1,2, which are not observable through neutrino oscillations. The probability that a neutrino created with a given flavor να is detected as a different flavor νβ after traveling a distance L in vacuum is given by [5, 6]
(2)Pνα→νβ=δαβ−4∑k>j3ℛeUαk∗UβkUαjUβj∗sin2Δmkj2L4E+2∑k>j3ℐmUαk∗UβkUαjUβj∗sinΔmkj2L2E,where E is the neutrino energy and Δmkj2≡mk2−mj2 are the differences of the squared masses of the definite-mass states k and j.In the period from the late 1990’s to early 2010’s, definitive experimental confirmation of neutrino oscillations was gathered from atmospheric [7, 8], solar (e.g., SNO [9, 10]), long baseline accelerator (K2K [11], MINOS [12], T2K [13]), and very long baseline reactor (Kam-LAND [14]) neutrino experiments. In 2012, the long baseline reactor neutrino experiments Daya Bay [15], RENO [16], and Double Chooz [17] reported the first measurement of the mixing angle θ13, by observing the disappearance of reactor antineutrinos ν¯e over distances of the order of 1 km, finding a nonzero value, and opening the door to studying CP violation in the neutrino sector. In later years, neutrino oscillation research has entered into a precision era, where the oscillation parameters can be determined with percent level precision from analyses that combine the results of many different experiments [18, 19]. The current focus of the field is primarily oriented to the determination of the CP-violating phase δCP, the neutrino mass ordering, and the octant of the angle θ23 (see for instance [20, 21] for recent experimental results and [22, 23] for current progress on new experimental efforts).In this work, we perform a combined analysis of the data from two specific data-taking periods of the Daya Bay and RENO experiments. For Daya Bay, we consider the 217 days of data, taken between December of 2011 and July of 2012 in the configuration with only 6 antineutrino detectors [24]. In the case of RENO, we consider the data from 500 live days, taken between August 2011 and January 2013 [25] with both the near and far detectors. These data sets have been chosen for the convenience and relative simplicity in reproducing their results from publicly available resources. Although more recent data are available, we have not considered them here. The aim of this work is to study the level of precision that can be attained by combining such older data sets, as well as to test our reproduction of the Daya Bay result with a full covariance matrix approach. In the following sections, we provide information about the experiments as well as a description of our analysis and results.


## 2. Antineutrinos from Nuclear Reactors

Typical commercial pressurized water reactors (PWRs) are copious sources ofν¯es, originating primarily in the beta decay of unstable fission products of the fissile isotopes 235U, 238U, 239Pu, and 241Pu, as well as from neutron capture in 238U. On average, each fission releases roughly 200 MeV of energy and produces 6 antineutrinos with energies below 10 MeV. Since the typical decay chain of the fission products has three consecutive beta decays, about 2×1020ν¯e per second per gigawatt of thermal power (GWth) [6] are isotropically emitted from the reactor core. The neutron capture contribution occurs at a smaller rate (0.6 per fission) and produces ν¯es with energies below 1.3 MeV [26].Reactor antineutrinos with energies>1.8MeV can be detected via the inverse beta decay (IBD) process
(3)ν¯e+p⟶e++n,by recording the delayed coincidence of the positron and neutron capture signals in, for example, a scintillating detector doped with a high neutron capture cross-section material, like gadolinium (Gd). The convolution of the reactor antineutrino flux and the IBD cross section gives an energy spectrum of the detected ν¯es with a peak around 3 MeV, and a cutoff at the 1.8 MeV threshold of the reaction.After the initial observation by RENO of a feature in theν¯e spectrum that has come to be known as the “5 MeV bump,” and its subsequent confirmation by Daya Bay, Double-Chooz, and other experiments, significant interest has arisen to try to explain it within the boundaries of nuclear physics, as well as through nonstandard particle physics (see [27–31] and references therein). The oscillation analyses developed by RENO and Daya Bay, which we reproduce here, assume that the bump is unrelated to the physics of neutrino oscillations and are mostly unaffected by this feature.


### 2.1. Brief Description of the Experiments

The Daya Bay experiment is located nearly 55 km northwest of Hong Kong; it uses the antineutrinos emitted by six functionally identical PWRs of 2.9 GWth each [15], two of them located in the Daya Bay Nuclear Power Plant (NPP) and four in the neighboring Ling Ao and Ling Ao-II NPPs. In Figure 1(a), the red circles represent the nuclear reactors, arranged in pairs (2 in Ling Ao, 2 in Ling Ao-II, and 2 in Daya Bay), and the blue circles represent the antineutrino detectors (ADs). In the data taking period used in this work, the detectors were distributed in three experimental halls (EH1, EH2, and EH3) as depicted in the figure. The three EHs are, respectively, under 250, 265, and 860 m.w.e. of overburden and are interconnected through internal tunnels, in order to shield the detectors from cosmic rays and other sources of radiation. The full 8 AD configuration was completed in 2012. Further details can be found in [32].Figure 1
Arrangement of the nuclear reactors, shown as red circles, and antineutrino detectors (ADs), shown as blue circles, for the Daya Bay (6 AD configuration) and RENO reactor neutrino experiments. Adapted from [33].
(a)
Daya Bay
(b)
RENOThe Reactor Experiment for Neutrino Oscillation (RENO) is located in the Hanbit (formerly Yonggwang) NPP in the southwest coast of South Korea, 250 km south of Seoul [34], and uses the ν¯e from six PWRs arranged in a line along the coast. The reactors produce a total of 16.4 GWth. RENO uses two detectors (near detector (ND) and far detector (FD)) to observe the produced ν¯es, as displayed in Figure 1(b), where the red circles represent the six reactors and the blue circles represent the two detectors. The ND (FD) is under 120 (450) m.w.e. of overburden [35]. The average distance from the reactors to the ND (FD) is 292 m (1380 m) [25].Daya Bay and RENO observeν¯e’s through the IBD reaction, Equation (3). Both experiments use a similar detector design with three concentric cylinders containing different liquids (see Figure 2). The inner-most cylinder is filled with a Gd-doped liquid scintillator (LS) and acts as the main target volume; the intermediate one, designed to efficiently detect gamma rays (gamma catcher), is filled with pure LS, and the outer one, whose inner walls are lined with photomultiplier tubes (PMTs), is filled with mineral oil, which acts as a buffer. The detectors are immersed in water pools whose walls are instrumented with PMTs and work as vetoes. Details of the detector design of each experiment can be found in [32] for DB and [36] for RENO.Figure 2
Schematic view of the Daya Bay ((a), adapted from [37]) and RENO ((b) [35]) antineutrino detectors, with the three concentric vessels (containing Gd-LS, LS, and MO) clearly identified. The interaction of a ν¯e with a proton in the target via IBD is shown inside the Daya Bay inner-most cylinder.
(a)
(b)In the target volume, a large amount of freely moving protonsp may interact with the antineutrinos coming from the reactors, producing positrons e+, which then annihilate with surrounding electrons generating two gamma rays (prompt signal). In addition, in the IBD process, a neutron n is also created; this thermalizes and is captured by a Gd nucleus, emitting more gamma rays (delayed signal). The time difference between these two signals is a few μs. A representation of the particle identification signal is shown inside the Daya Bay detector in the top panel of Figure 2. Once the detected signals are collected, specialized selection criteria are applied by the experiments to estimate the observed number of ν¯e events and background rates in each detector (for detailed information about this process see, for instance, [15, 16, 38]).


### 2.2. Input to Our Studies

The prompt reconstructed energy,Ep, distributions used in the analyses are shown in Figure 3. For Daya Bay, we digitized the data and no-oscillation distributions from Figure 2 in Ref. [24] and assume that all detectors in the same EH have the same distribution. Note that the predicted Daya Bay distributions already account for the 5 MeV bump. For RENO, we digitized the data and best-fit distributions from Figure 26 in Ref. [25]. The no-oscillation distributions in the near and far RENO detectors were constructed by removing from the best-fit spectrum, bin by bin, the effect of the oscillations with the help of a sample of simulated neutrino events (see Section 2.3). Note that in the RENO case, the predicted distributions do not include the 5 MeV bump; however, their spectral analysis, based on a far-to-near ratio, described later, will prove to be insensitive to this effect.Figure 3
Background-subtracted reconstructed positron prompt energy distributions for Daya Bay (a) and RENO (b). Black crosses are data; blue (red) line histograms show the MC expected no-oscillation (best fit) spectra. The background histograms are shown in magenta. See text for details. For Daya Bay, the distributions are absolutely normalized and account for the 5 MeV bump effect. For RENO, the predicted ND and FD distributions are normalized such that the best fit and data histograms have equal areas outside the 5 MeV bump region (3.8MeV<Ep<6.4MeV). The bottom-right plot (RENO) compares the observed spectrum in the FD with no oscillation and the best-fit predictions obtained from the measurement at the ND, which account for the effect of the 5 MeV bump.
(a)
(b)In order to normalize the event rates and energy distributions, we collected the information from tables found in Refs. [24, 25], which we have summarized here in Tables 1 and 2, for Daya Bay and RENO, respectively. In these tables, we have included the estimates of the total IBD rates without oscillations used in our simulation for each experiment. In our RENO simulation, we set the total predicted IBD rate at the best fit to the observed value and used the approximation of a common detection efficiency for the near and far detectors. In addition, for the Daya Bay spectral analysis, we digitized the full systematic error correlation matrix from Ref. [39], the total systematic errors from Figure 2 in Ref. [40], and used this information to construct the full covariance matrix, as will be described in Section 3 below.Table 1
Daya Bay IBD candidates and expected IBD rates and backgrounds per day in the 217-day sample used in this work. The DAQ live time and the muon veto and multiplicity-cut efficiencies (∈μ⋅∈m) are also reported for each AD. The last line shows the expected IBD rate without oscillations. Adapted from [13].
Experimental hall
EH1
EH2
EH3
Detector
AD1
AD2
AD3
AD4
AD5
AD6
IBD candidates/day
530.21±1.67
536.75±1.68
489.93±1.61
73.58±0.62
73.21±0.62
72.35±0.62
DAQ live time (days)
191.001
189.645
189.779
∈μ⋅∈m
0.7957
0.7927
0.8282
0.9577
0.95689
0.9566
Total background/day
13.20±0.98
13.01±0.98
9.57±0.71
3.52±0.14
3.48±0.14
3.43±0.14
IBD rates/day (best fit)
653.30±2.31
664.15±2.33
581.97±2.07
73.31±0.66
73.03±0.66
72.20±0.66
No-osc. IBD rates/day
665.11
675.95
593.99
78.94
78.64
77.74Table 2
RENO IBD rates and backgrounds per day in the 500 live days used in this work. DAQ live times are also reported for the far and near detectors. A common detection efficiency was used in our simulation. Adapted from [41].
Detector
Near
Far
IBD rate (background sub.)
616.67±1.44
61.24±0.42
DAQ live time (days)
458.49
489.93
Detection efficiency
0.7644
0.7644
Total background
17.54±0.83
3.14±0.23
No-osc. IBD rate
623.26
65.57


### 2.3. Simulation of Neutrino Events

We simulate neutrino events traveling the different baselines available between the various reactors and detectors in each experiment by constructing the probability that a neutrino leaves a particular reactorr and arrives at a specific detector d. For RENO (6 reactors and 2 detectors), there are 12 different baselines, while for Daya Bay (6 detectors and 6 reactors), there are 36 different baselines, in the 6 AD configuration considered here. This probability is calculated as follows:
(4)wrd=NPrthMd4πLrd2,where Prth is the thermal power of reactor r, Md is the mass of the fiducial volume of detector d, and Lrd is the baseline distance between reactor r and detector d; N is a normalization constant making the sum ∑rdwrd=1.The baseline lengths and detector fiducial volume masses were extracted from Ref. [32] for the Daya Bay analysis, and from Ref. [34] for the RENO analysis. Figure 4 shows the probability distributions for a neutrino to travel along each available baseline. Note that each baseline index (1–12 for RENO and 1-36 for Daya Bay) uniquely identifies a reactor-detector pair. The histograms give the probability that a neutrino in the experiment was produced in a particular reactor and observed in a particular detector. As expected, near detectors have larger probabilities than far detectors, since the closer the detector is to a reactor, the more antineutrinos are detected.Figure 4
Baseline probability distributions for Daya Bay (a) and RENO (b), obtained from the simulation of 107 neutrino events. Histograms are normalized to unit area.
(a)
(b)The neutrino energy for each simulated event is then assigned using the well-known relation [13]
(5)Eν=Ep+E¯n+0.78MeV,where E¯n is the average energy taken by the neutron (∼10 keV). A small multiplicative correction factor f=1.03 was applied to EpEp⟶Ep×f in our RENO simulation to better reproduce the results for this experiment. Since we do not have access to the no-oscillated RENO spectra, this factor, which is degenerate with Δmee2, serves the purpose to make our overly simplistic simulation resemble more closely that developed by the collaboration. While this neutrino energy is rather a reconstructed quantity and not the actual true energy of the event, the energy resolution effect will be neglected, as it is reported to be smaller than the bin size.As mentioned in Section2.2, Refs. [24, 25] only report the number of either expected or observed IBD events in each detector, with the effect of oscillations at the best-fit. However, knowledge of the number of IBD candidate events expected without oscillations is required. We estimated the number of events without oscillations dividing the number of oscillated events in a given bin, or full spectrum, by the average best-fit oscillation probability of all the events in said bin or spectrum. Average oscillation probabilities were calculated by applying the oscillation probability in Equation (7) to a sample of 107 simulated neutrino IBD events whose true energy Eν, prompt positron energy Ep, and baseline L are assigned as follows: first, a baseline L is randomly sampled from the distribution of baselines in Figure 4; this uniquely determines the reactor-detector pair for the event. A random Gaussian fluctuation (σ=1m) is added to approximately incorporate the reactor and detector sizes. The prompt positron energy is then sampled from the Ep distribution corresponding to the chosen detector (and corrected by multiplying it by f). This Ep value is then used to calculate the true neutrino energy Eν using Equation (5).


## 2.1. Brief Description of the Experiments

The Daya Bay experiment is located nearly 55 km northwest of Hong Kong; it uses the antineutrinos emitted by six functionally identical PWRs of 2.9 GWth each [15], two of them located in the Daya Bay Nuclear Power Plant (NPP) and four in the neighboring Ling Ao and Ling Ao-II NPPs. In Figure 1(a), the red circles represent the nuclear reactors, arranged in pairs (2 in Ling Ao, 2 in Ling Ao-II, and 2 in Daya Bay), and the blue circles represent the antineutrino detectors (ADs). In the data taking period used in this work, the detectors were distributed in three experimental halls (EH1, EH2, and EH3) as depicted in the figure. The three EHs are, respectively, under 250, 265, and 860 m.w.e. of overburden and are interconnected through internal tunnels, in order to shield the detectors from cosmic rays and other sources of radiation. The full 8 AD configuration was completed in 2012. Further details can be found in [32].Figure 1
Arrangement of the nuclear reactors, shown as red circles, and antineutrino detectors (ADs), shown as blue circles, for the Daya Bay (6 AD configuration) and RENO reactor neutrino experiments. Adapted from [33].
(a)
Daya Bay
(b)
RENOThe Reactor Experiment for Neutrino Oscillation (RENO) is located in the Hanbit (formerly Yonggwang) NPP in the southwest coast of South Korea, 250 km south of Seoul [34], and uses the ν¯e from six PWRs arranged in a line along the coast. The reactors produce a total of 16.4 GWth. RENO uses two detectors (near detector (ND) and far detector (FD)) to observe the produced ν¯es, as displayed in Figure 1(b), where the red circles represent the six reactors and the blue circles represent the two detectors. The ND (FD) is under 120 (450) m.w.e. of overburden [35]. The average distance from the reactors to the ND (FD) is 292 m (1380 m) [25].Daya Bay and RENO observeν¯e’s through the IBD reaction, Equation (3). Both experiments use a similar detector design with three concentric cylinders containing different liquids (see Figure 2). The inner-most cylinder is filled with a Gd-doped liquid scintillator (LS) and acts as the main target volume; the intermediate one, designed to efficiently detect gamma rays (gamma catcher), is filled with pure LS, and the outer one, whose inner walls are lined with photomultiplier tubes (PMTs), is filled with mineral oil, which acts as a buffer. The detectors are immersed in water pools whose walls are instrumented with PMTs and work as vetoes. Details of the detector design of each experiment can be found in [32] for DB and [36] for RENO.Figure 2
Schematic view of the Daya Bay ((a), adapted from [37]) and RENO ((b) [35]) antineutrino detectors, with the three concentric vessels (containing Gd-LS, LS, and MO) clearly identified. The interaction of a ν¯e with a proton in the target via IBD is shown inside the Daya Bay inner-most cylinder.
(a)
(b)In the target volume, a large amount of freely moving protonsp may interact with the antineutrinos coming from the reactors, producing positrons e+, which then annihilate with surrounding electrons generating two gamma rays (prompt signal). In addition, in the IBD process, a neutron n is also created; this thermalizes and is captured by a Gd nucleus, emitting more gamma rays (delayed signal). The time difference between these two signals is a few μs. A representation of the particle identification signal is shown inside the Daya Bay detector in the top panel of Figure 2. Once the detected signals are collected, specialized selection criteria are applied by the experiments to estimate the observed number of ν¯e events and background rates in each detector (for detailed information about this process see, for instance, [15, 16, 38]).


## 2.2. Input to Our Studies

The prompt reconstructed energy,Ep, distributions used in the analyses are shown in Figure 3. For Daya Bay, we digitized the data and no-oscillation distributions from Figure 2 in Ref. [24] and assume that all detectors in the same EH have the same distribution. Note that the predicted Daya Bay distributions already account for the 5 MeV bump. For RENO, we digitized the data and best-fit distributions from Figure 26 in Ref. [25]. The no-oscillation distributions in the near and far RENO detectors were constructed by removing from the best-fit spectrum, bin by bin, the effect of the oscillations with the help of a sample of simulated neutrino events (see Section 2.3). Note that in the RENO case, the predicted distributions do not include the 5 MeV bump; however, their spectral analysis, based on a far-to-near ratio, described later, will prove to be insensitive to this effect.Figure 3
Background-subtracted reconstructed positron prompt energy distributions for Daya Bay (a) and RENO (b). Black crosses are data; blue (red) line histograms show the MC expected no-oscillation (best fit) spectra. The background histograms are shown in magenta. See text for details. For Daya Bay, the distributions are absolutely normalized and account for the 5 MeV bump effect. For RENO, the predicted ND and FD distributions are normalized such that the best fit and data histograms have equal areas outside the 5 MeV bump region (3.8MeV<Ep<6.4MeV). The bottom-right plot (RENO) compares the observed spectrum in the FD with no oscillation and the best-fit predictions obtained from the measurement at the ND, which account for the effect of the 5 MeV bump.
(a)
(b)In order to normalize the event rates and energy distributions, we collected the information from tables found in Refs. [24, 25], which we have summarized here in Tables 1 and 2, for Daya Bay and RENO, respectively. In these tables, we have included the estimates of the total IBD rates without oscillations used in our simulation for each experiment. In our RENO simulation, we set the total predicted IBD rate at the best fit to the observed value and used the approximation of a common detection efficiency for the near and far detectors. In addition, for the Daya Bay spectral analysis, we digitized the full systematic error correlation matrix from Ref. [39], the total systematic errors from Figure 2 in Ref. [40], and used this information to construct the full covariance matrix, as will be described in Section 3 below.Table 1
Daya Bay IBD candidates and expected IBD rates and backgrounds per day in the 217-day sample used in this work. The DAQ live time and the muon veto and multiplicity-cut efficiencies (∈μ⋅∈m) are also reported for each AD. The last line shows the expected IBD rate without oscillations. Adapted from [13].
Experimental hall
EH1
EH2
EH3
Detector
AD1
AD2
AD3
AD4
AD5
AD6
IBD candidates/day
530.21±1.67
536.75±1.68
489.93±1.61
73.58±0.62
73.21±0.62
72.35±0.62
DAQ live time (days)
191.001
189.645
189.779
∈μ⋅∈m
0.7957
0.7927
0.8282
0.9577
0.95689
0.9566
Total background/day
13.20±0.98
13.01±0.98
9.57±0.71
3.52±0.14
3.48±0.14
3.43±0.14
IBD rates/day (best fit)
653.30±2.31
664.15±2.33
581.97±2.07
73.31±0.66
73.03±0.66
72.20±0.66
No-osc. IBD rates/day
665.11
675.95
593.99
78.94
78.64
77.74Table 2
RENO IBD rates and backgrounds per day in the 500 live days used in this work. DAQ live times are also reported for the far and near detectors. A common detection efficiency was used in our simulation. Adapted from [41].
Detector
Near
Far
IBD rate (background sub.)
616.67±1.44
61.24±0.42
DAQ live time (days)
458.49
489.93
Detection efficiency
0.7644
0.7644
Total background
17.54±0.83
3.14±0.23
No-osc. IBD rate
623.26
65.57


## 2.3. Simulation of Neutrino Events

We simulate neutrino events traveling the different baselines available between the various reactors and detectors in each experiment by constructing the probability that a neutrino leaves a particular reactorr and arrives at a specific detector d. For RENO (6 reactors and 2 detectors), there are 12 different baselines, while for Daya Bay (6 detectors and 6 reactors), there are 36 different baselines, in the 6 AD configuration considered here. This probability is calculated as follows:
(4)wrd=NPrthMd4πLrd2,where Prth is the thermal power of reactor r, Md is the mass of the fiducial volume of detector d, and Lrd is the baseline distance between reactor r and detector d; N is a normalization constant making the sum ∑rdwrd=1.The baseline lengths and detector fiducial volume masses were extracted from Ref. [32] for the Daya Bay analysis, and from Ref. [34] for the RENO analysis. Figure 4 shows the probability distributions for a neutrino to travel along each available baseline. Note that each baseline index (1–12 for RENO and 1-36 for Daya Bay) uniquely identifies a reactor-detector pair. The histograms give the probability that a neutrino in the experiment was produced in a particular reactor and observed in a particular detector. As expected, near detectors have larger probabilities than far detectors, since the closer the detector is to a reactor, the more antineutrinos are detected.Figure 4
Baseline probability distributions for Daya Bay (a) and RENO (b), obtained from the simulation of 107 neutrino events. Histograms are normalized to unit area.
(a)
(b)The neutrino energy for each simulated event is then assigned using the well-known relation [13]
(5)Eν=Ep+E¯n+0.78MeV,where E¯n is the average energy taken by the neutron (∼10 keV). A small multiplicative correction factor f=1.03 was applied to EpEp⟶Ep×f in our RENO simulation to better reproduce the results for this experiment. Since we do not have access to the no-oscillated RENO spectra, this factor, which is degenerate with Δmee2, serves the purpose to make our overly simplistic simulation resemble more closely that developed by the collaboration. While this neutrino energy is rather a reconstructed quantity and not the actual true energy of the event, the energy resolution effect will be neglected, as it is reported to be smaller than the bin size.As mentioned in Section2.2, Refs. [24, 25] only report the number of either expected or observed IBD events in each detector, with the effect of oscillations at the best-fit. However, knowledge of the number of IBD candidate events expected without oscillations is required. We estimated the number of events without oscillations dividing the number of oscillated events in a given bin, or full spectrum, by the average best-fit oscillation probability of all the events in said bin or spectrum. Average oscillation probabilities were calculated by applying the oscillation probability in Equation (7) to a sample of 107 simulated neutrino IBD events whose true energy Eν, prompt positron energy Ep, and baseline L are assigned as follows: first, a baseline L is randomly sampled from the distribution of baselines in Figure 4; this uniquely determines the reactor-detector pair for the event. A random Gaussian fluctuation (σ=1m) is added to approximately incorporate the reactor and detector sizes. The prompt positron energy is then sampled from the Ep distribution corresponding to the chosen detector (and corrected by multiplying it by f). This Ep value is then used to calculate the true neutrino energy Eν using Equation (5).


## 3. Oscillation Analysis

Neutrino oscillations in long baseline reactor neutrino experiments manifest themselves as the disappearance of electron antineutrinos with energies between 2 and 6 MeV over distances of the order of 1 km. In this case, the survival probabilityPν¯e→ν¯e in Equation (2) is well approximated by
(6)Pν¯e→ν¯e=1−cos4θ13sin22θ12sin2Δ21−sin22θ13cos2θ12sin2Δ31+sin2θ12sin2Δ32,with Δij=1.267Δmij2L/Eν, and Δmij2 is given in eV2, L (in m) is the distance between the reactor and the detector, and Eν (in MeV) is the neutrino energy. Given that Δm212≪Δm312≈Δm322, the ν¯e oscillation is mainly driven by Δ31, and Equation (6) naturally leads to the definition of an effective squared-mass difference Δmee2 such that sin2Δee=cos2θ12sin2Δ31+sin2θ12sin2Δ32 [24], so that the survival probability in Equation (6) can be reduced to
(7)Pν¯e→ν¯e=1−cos4θ13sin22θ12sin2Δ21−sin22θ13sin2Δee.For the purpose of reproducing the Daya Bay results, in what follows, we have usedsin22θ12=0.857, Δm212=7.50×10−5eV2, for the spectral analysis, and also Δm322=2.32×10−3eV2, for the rate-only analysis, as in Ref. [24]. For the reproduction of the RENO results, we used sin22θ12=0.846, and Δm212=7.53×10−5eV2, for the spectral analysis, and in addition Δm322=2.49×10−3eV2 for the rate-only analysis, as in Ref. [25].It has been pointed out that the definition ofΔmee2 used here, besides being L/E dependent, is discontinuous at 0.5 km/MeV [42], and better definition can be considered, such as the weighted average of Δm312 and Δm322. In the interest of attempting to reproduce the original results by Daya Bay and RENO, we will keep the definition in Ref. [24], as this has no critical impact for our analysis.We now present the procedure to estimate the oscillation parameters which best fit the data, preforming two types of approaches: a rate-only analysis and a spectral analysis.


### 3.1. Rate-Only Analysis

Using only the information of the total event rates reported in Tables1 and 2, and the measured values of the oscillation parameters sin22θ12, Δm212, and Δm322, we extract the value of sin22θ13. Since the shape of the spectrum is not considered, no information about Δmee2 is obtained from this analysis.For Daya Bay, we follow Ref. [15] and define our χ2 statistic as
(8)χ2=∑d=16Md−Td1+ε+∑rωrdαr+εd+ηd2Md+Bd+∑rαr2σr2+∑d=16εd2σd2+ηd2σB2.Here,Md is the number of IBD-observed events in the dth detector after background subtraction, Bd is the background rate, and Td is the predicted number of events considering neutrino oscillations through Equation (7), in the total DAQ live time (Table 1). The quantity wrd is the fraction of events produced in reactor r which contribute to detector d, considering the travel distance and the neutrino flux, which corresponds precisely with the baseline probability of Figure 4. The χ2 in (8) is penalized by the inclusion of 18 pull terms αr, εd, and ηd (with r,d=1,⋯,6), characterizing the systematic errors affecting the measurement. As reported by the collaboration in Ref. [15], σr (0.8%) and σd (0.2%) are the uncorrelated reactor and detector uncertainties, respectively, and σB is the corresponding background uncertainty. The additional parameter ε is included as a normalization factor which accounts for possible differences between the observation and the prediction, and it is included as a free parameter.The interval0<sin22θ13<0.2 is split in 200 uniform steps. For each point, a full minimization over the 18 pull terms and the parameter ε is performed to obtain the value of the marginalized χ2 statistic. Minimization of the marginalized χ2 gives the best-fit value sin22θ13=0.090−0.009+0.010 at 1σ C.L.A similar rate-only analysis was performed to the RENO data, considering the three-neutrino oscillation model described by Equation (7). In this case, we follow Ref. [16] and define the χ2 statistic as
(9)χ2=∑d=N,FNobsd+bd−1+a+ξd∑r1+frNexpd,r2Nobs+∑d=N,Fξd2σdξ2+bd2σdb2+∑r=16frσr2,where Nobsd is the number of observed events after background subtraction for the near (N) and far (F) detectors in the total DAQ live time (Table 2); Nexpd,r is the number of expected events in detector d coming from reactor r, including the detection efficiency and the effect of oscillations. Here, a is a freely varying normalization factor, σbd, from Table 2, are the background uncertainties associated with the pull term parameters bd, and σr (0.9%) and σdξ (0.2%) are the uncorrelated reactor and detector uncertainties, associated with the pull term parameters fr and ξd, respectively. Using a similar minimization procedure to the one used for the Daya Bay analysis, marginalizing over a and the pull term parameters, we obtain sin22θ13=0.088−0.013+0.010 at 1σ as the value which best fit to the RENO data.


### 3.2. Spectral Analysis

Here, besides the normalization information, the shape of energy distribution of the observed IBD events will be used to extract the value ofΔmee2 along with sin22θ13. The best fit is found by minimizing a suitable χ2 statistic defined over a uniformly spaced 100×100 grid in the sin22θ13vs.Δmee2 space.


#### 3.2.1. Daya Bay

For Daya Bay, we follow Ref. [43] and define
(10)χ2=∑i,jNiobs−NiexpTVij−1Njobs−Njexp,where the indices i and j run over 156 bins corresponding to the concatenation of the prompt energy distributions of the 6 ADs, with 26 bins each. NiobsNiexp is the number of observed (expected) events in the ith energy bin, and Vij are the elements of the total covariance matrix expressed in the same binning. Niexp depends on the oscillation parameters sin22θ13 and Δmee2. The covariance matrix, V=Vstat+Vsyst, contains all sources of statistic and systematic errors affecting the experiment. The systematic error component was calculated from the full correlation matrix reported in [39], including the signal, background, and reactor core errors, and the total systematic uncertainties presented in Figure 2 of Ref. [40]. This assumption proved to be a reasonable approximation to the total systematic errors and correlations in the data set used for this analysis.The systematic error correlation matrix in [39] is a 222×222 matrix, since the energy spectrum in each AD used therein has 37 (nonuniform) bins in the energy range 0.7-12 MeV. However, the Daya Bay data set considered here used energy distributions with 26 bins, in the same energy range, hence having different bin boundaries (see Figure 5). In order to cast this correlation matrix in the form of a 156×156 matrix, we implemented a rebinning procedure based on the diagonal blocks (37×37 bins per AD) of the original matrix in which we sampled energy distributions consistent with the original matrix [44], and for each one, sampled 106 energy values which were then filled into a rebinned histogram (26×26 bins per AD). The resulting 1000 rebinned distributions were used to recalculate the correlation matrix with the desired binning. The full 6 AD correlation matrix, ρsyst, was constructed assuming a correlation among the 6 detectors encoded in a 6×6 matrix (see middle plot in Figure 6) which was adjusted so that the overall features of the original matrix could be reproduced. Our rebinned correlation matrix is shown in the top panel of Figure 6. The 156×156 elements of the total systematic covariance matrix Vsyst were then computed as
(11)Vijsyst=ρijsyst×Niexpσimod26×Njexpσjmod26(no summation over repeated indices) where the σk, k=1,⋯,26, are the total fractional systematic uncertainties, shown in the bottom panel of Figure 6, extrapolated from [40] up to 12 MeV. Despite being a rough approximation to the true error matrix used by the collaboration, as we will see, our results agree reasonably well with those reported by Daya Bay.Figure 5
Energy range and bin boundaries for the Daya Bay analysis in [39] (black lines and numbers), and for the analysis presented in this work (blue lines and numbers).Figure 6
(a) Full systematic error (signal, background, and reactor cores) correlation matrix used in the Daya Bay analysis with 26 bins per AD (signal, backgrounds, and reactor cores). Each square block represents the correlation among the bins of two ADs. (b) Correlation matrix for the 6 ADs. (c) Total fractional systematic uncertainties in the prompt energy distribution, assumed equal for all ADs.
(a)
(b)
(c)We found that the oscillation parameters which best reproduce the data aresin22θ13=0.091−0.009+0.012 and Δmee2=2.60−0.22+0.18×10−3eV2 at 1σ C.L., with χmin2/NDF=59.91/154. This result is shown in Figure 7, together with the 68.27%, 95.45%, and 99.73% C.L. allowed regions for the oscillation parameters space. In the upper (right) panel of Figure 7, we show the Δχ2 marginalization over Δmee2sin22θ13, where the minimum of the curve Δχ2=0 points to the best-fit value. We have included here the result of the rate-only analysis (dash-dotted line in the top panel) for comparison purposes. In the marginalization plots, the horizontal (vertical) lines indicate the one-dimensional allowed regions for the two parameters at the same C.L. as the 2D plot.Figure 7
Allowed regions in thesin22θ13,Δmee2 space at a 68.27%, 95.45%, and 99.73% C.L. from our Daya Bay analysis. The best fit is marked by a “+” in the central plot, where the published Daya Bay contours [15] are also included. The marginalized Δχ2 over the oscillation parameters are also shown in the upper and right plots, including the rate-only result (dot-dashed line) for the mixing angle, and the 99.99% C.L. line (solid dark blue).Although the best fit is very well recovered, our contours are slightly wider than the published ones (dotted line in Figure7) towards the higher sin22θ13 values, and shorter towards the lower Δmee2. Despite our efforts to reproduce the full covariance matrix for this measurement, several manipulations had to be implemented in order to rebin the matrix and guarantee its positive definiteness, which may have introduced distortions. Nonetheless, we consider that our result captures the main features of the analysis and gives a good approximation to the confidence regions for the parameters. A cross-check calculation using a χ2 with pull terms produced contours with similar characteristics.


#### 3.2.2. RENO

Following [25], we define the χ2 statistic for the RENO spectral analysis as
(12)χ2=∑i=1NbinsOiF/N−TiF/N2UiF/N+∑d=N,Fbdσbkdd2+eσscale2+∈σeff2.In this case,OiF/N is the ratio of the observed IBD candidate events at the far detector over those observed at the near detector for the ith energy bin; TiF/N=TiF/Nbd,e,∈,θ13,Δmee2 is the corresponding ratio of expected events, and UiF/N is the statistical uncertainty associated with OiF/N.The best fit found after minimizing theχ2 statistic with respect to the oscillation parameters, marginalizing over bd, e, and ∈ (the pull term parameters), is sin22θ13=0.083−0.012+0.010 and Δmee2=2.64−0.27+0.21×10−3eV2 at 1σ C.L., with χmin2/NDF=20.92/25. The results of this analysis are presented in Figure 8, where the allowed regions in the sin22θ13,Δmee2 parameter space are shown, together with the best-fit point. As in the case for the Daya Bay analysis, the 1-dimensional marginalized distributions are also shown for each of the oscillation parameters, and the result obtained from the rate-only analysis (dashed line in top panel) for sin22θ13 is included for comparison. We have also included the contours and best fit obtained by RENO [25] in Figure 8 (dotted line) which show good agreement with our analysis.Figure 8
Allowed regions in the (sin22θ13,Δmee2) space at a 68.27%, 95.45%, and 99.73% C.L. from our RENO analysis. The best fit is marked by a “+” in the central plot, where the published RENO contours [25] are also included. The marginalized Δχ2 over the oscillation parameters are also shown in the upper and right plots, including the rate-only result (dot-dashed line) for the mixing angle, and the 99.99% C.L. line (solid dark blue).The bottom-right plot in Figure3 compares the FD data to the no oscillations and best-fit predictions obtained from the measured spectrum at the ND. The agreement between the data and the best-fit prediction is very good and demonstrates that the near-to-far ratio technique used in the RENO analysis is insensitive to the presence of the 5 MeV bump.


### 3.3. Combined Analysis

Finally, we performed a combined analysis of the two data sets by considering theχ2 statistic
(13)χcomb2=ΔχDayaBay2+ΔχRENO2.This definition will answer the specific question: how probable is it that both experimental results come from the same underlying oscillation model? [45]. Both the rate-only and the spectral analysis were performed using the corresponding Δχ2=χ2−χmin2 statistic appropriate for each case.For the combined rate-only analysis, we found that the data is best described with a value ofsin22θ13=0.088−0.006+0.008. This result is shown in Figure 9, where Δχ2=χcomb2−χcombmin2 is plotted as a function of sin22θ13 (solid purple line). We also plot here the rate-only results from the independent analyses of Daya Bay (dotted line) and RENO (dash-dotted line). Horizontal colored lines are drawn to mark the allowed intervals for the parameter at 1-4σ C.L., which are smaller for the combined analysis, indicating the expected enhancement in significance.Figure 9
ObtainedΔχ2 distribution for sin22θ13 from the Daya Bay+RENO rate-only combined analysis (full purple line), together with the results from Daya Bay (dotted line) and RENO (dash-dotted line), separately.For the spectral combined analysis, theχcomb2 is also built as the sum of the corresponding statistics used for Daya Bay and RENO, as in (13). The minimization of such a function gives the values of the oscillation parameters which produce the best fit to the data: sin22θ13=0.087−0.008+0.007, Δmee2=2.59−0.14+0.15×10−3eV21σ, with χcombmin2/NDF=0.18/2. Together with the best-fit point, Figure 10 shows the 68.27%, 95.45%, and 99.73% C.L. allowed regions in the studied parameter space.Figure 10
2D regions (1, 2, and 3σ) and 1D intervals (1, 2, 3, and 4σ) allowed for the oscillation parameters obtained from the combination of Daya Bay and RENO data, using the spectral analysis. The best fit is marked by a “+” in the central plot. The results from the independent analyses (RENO dot-dashed; Daya Bay dotted) are shown for comparison.We also show in Figure10 (top and right panels) the 1-dimensional Δχ2 distributions for each oscillation parameter, marginalizing over the other one, where we have included the results obtained separately for Daya Bay (dotted line) and RENO (dash-dotted line). Clearly, the combined result (solid line) produces smaller intervals for sin22θ13 and Δmee2.Finally, using the prescription described in Ref. [45], we evaluate the compatibility between the two data sets, by calculating the parameter goodness (PG) as the χ-squared probability Prob(χcombmin2; Pc), where Pc=2 is the number of parameters coupling the two data sets. In this case, we obtain a compatibility of PG=91.5%. We note that the small discrepancies between our results for Daya Bay and the published ones may be leading us to a different compatibility level from what could be obtained with the official analysis results. Our wider contour in sin22θ13 tends to increase the compatibility, while being smaller along Δmee2 tends to reduce it. However, with regard to the question at the beginning of this section, the compatibility found here allows us to state that both the Daya Bay and RENO data considered in this work are well described by the same neutrino oscillation model, represented by (7), with the oscillation parameters found in the combined analysis.


## 3.1. Rate-Only Analysis

Using only the information of the total event rates reported in Tables1 and 2, and the measured values of the oscillation parameters sin22θ12, Δm212, and Δm322, we extract the value of sin22θ13. Since the shape of the spectrum is not considered, no information about Δmee2 is obtained from this analysis.For Daya Bay, we follow Ref. [15] and define our χ2 statistic as
(8)χ2=∑d=16Md−Td1+ε+∑rωrdαr+εd+ηd2Md+Bd+∑rαr2σr2+∑d=16εd2σd2+ηd2σB2.Here,Md is the number of IBD-observed events in the dth detector after background subtraction, Bd is the background rate, and Td is the predicted number of events considering neutrino oscillations through Equation (7), in the total DAQ live time (Table 1). The quantity wrd is the fraction of events produced in reactor r which contribute to detector d, considering the travel distance and the neutrino flux, which corresponds precisely with the baseline probability of Figure 4. The χ2 in (8) is penalized by the inclusion of 18 pull terms αr, εd, and ηd (with r,d=1,⋯,6), characterizing the systematic errors affecting the measurement. As reported by the collaboration in Ref. [15], σr (0.8%) and σd (0.2%) are the uncorrelated reactor and detector uncertainties, respectively, and σB is the corresponding background uncertainty. The additional parameter ε is included as a normalization factor which accounts for possible differences between the observation and the prediction, and it is included as a free parameter.The interval0<sin22θ13<0.2 is split in 200 uniform steps. For each point, a full minimization over the 18 pull terms and the parameter ε is performed to obtain the value of the marginalized χ2 statistic. Minimization of the marginalized χ2 gives the best-fit value sin22θ13=0.090−0.009+0.010 at 1σ C.L.A similar rate-only analysis was performed to the RENO data, considering the three-neutrino oscillation model described by Equation (7). In this case, we follow Ref. [16] and define the χ2 statistic as
(9)χ2=∑d=N,FNobsd+bd−1+a+ξd∑r1+frNexpd,r2Nobs+∑d=N,Fξd2σdξ2+bd2σdb2+∑r=16frσr2,where Nobsd is the number of observed events after background subtraction for the near (N) and far (F) detectors in the total DAQ live time (Table 2); Nexpd,r is the number of expected events in detector d coming from reactor r, including the detection efficiency and the effect of oscillations. Here, a is a freely varying normalization factor, σbd, from Table 2, are the background uncertainties associated with the pull term parameters bd, and σr (0.9%) and σdξ (0.2%) are the uncorrelated reactor and detector uncertainties, associated with the pull term parameters fr and ξd, respectively. Using a similar minimization procedure to the one used for the Daya Bay analysis, marginalizing over a and the pull term parameters, we obtain sin22θ13=0.088−0.013+0.010 at 1σ as the value which best fit to the RENO data.


## 3.2. Spectral Analysis

Here, besides the normalization information, the shape of energy distribution of the observed IBD events will be used to extract the value ofΔmee2 along with sin22θ13. The best fit is found by minimizing a suitable χ2 statistic defined over a uniformly spaced 100×100 grid in the sin22θ13vs.Δmee2 space.


### 3.2.1. Daya Bay

For Daya Bay, we follow Ref. [43] and define
(10)χ2=∑i,jNiobs−NiexpTVij−1Njobs−Njexp,where the indices i and j run over 156 bins corresponding to the concatenation of the prompt energy distributions of the 6 ADs, with 26 bins each. NiobsNiexp is the number of observed (expected) events in the ith energy bin, and Vij are the elements of the total covariance matrix expressed in the same binning. Niexp depends on the oscillation parameters sin22θ13 and Δmee2. The covariance matrix, V=Vstat+Vsyst, contains all sources of statistic and systematic errors affecting the experiment. The systematic error component was calculated from the full correlation matrix reported in [39], including the signal, background, and reactor core errors, and the total systematic uncertainties presented in Figure 2 of Ref. [40]. This assumption proved to be a reasonable approximation to the total systematic errors and correlations in the data set used for this analysis.The systematic error correlation matrix in [39] is a 222×222 matrix, since the energy spectrum in each AD used therein has 37 (nonuniform) bins in the energy range 0.7-12 MeV. However, the Daya Bay data set considered here used energy distributions with 26 bins, in the same energy range, hence having different bin boundaries (see Figure 5). In order to cast this correlation matrix in the form of a 156×156 matrix, we implemented a rebinning procedure based on the diagonal blocks (37×37 bins per AD) of the original matrix in which we sampled energy distributions consistent with the original matrix [44], and for each one, sampled 106 energy values which were then filled into a rebinned histogram (26×26 bins per AD). The resulting 1000 rebinned distributions were used to recalculate the correlation matrix with the desired binning. The full 6 AD correlation matrix, ρsyst, was constructed assuming a correlation among the 6 detectors encoded in a 6×6 matrix (see middle plot in Figure 6) which was adjusted so that the overall features of the original matrix could be reproduced. Our rebinned correlation matrix is shown in the top panel of Figure 6. The 156×156 elements of the total systematic covariance matrix Vsyst were then computed as
(11)Vijsyst=ρijsyst×Niexpσimod26×Njexpσjmod26(no summation over repeated indices) where the σk, k=1,⋯,26, are the total fractional systematic uncertainties, shown in the bottom panel of Figure 6, extrapolated from [40] up to 12 MeV. Despite being a rough approximation to the true error matrix used by the collaboration, as we will see, our results agree reasonably well with those reported by Daya Bay.Figure 5
Energy range and bin boundaries for the Daya Bay analysis in [39] (black lines and numbers), and for the analysis presented in this work (blue lines and numbers).Figure 6
(a) Full systematic error (signal, background, and reactor cores) correlation matrix used in the Daya Bay analysis with 26 bins per AD (signal, backgrounds, and reactor cores). Each square block represents the correlation among the bins of two ADs. (b) Correlation matrix for the 6 ADs. (c) Total fractional systematic uncertainties in the prompt energy distribution, assumed equal for all ADs.
(a)
(b)
(c)We found that the oscillation parameters which best reproduce the data aresin22θ13=0.091−0.009+0.012 and Δmee2=2.60−0.22+0.18×10−3eV2 at 1σ C.L., with χmin2/NDF=59.91/154. This result is shown in Figure 7, together with the 68.27%, 95.45%, and 99.73% C.L. allowed regions for the oscillation parameters space. In the upper (right) panel of Figure 7, we show the Δχ2 marginalization over Δmee2sin22θ13, where the minimum of the curve Δχ2=0 points to the best-fit value. We have included here the result of the rate-only analysis (dash-dotted line in the top panel) for comparison purposes. In the marginalization plots, the horizontal (vertical) lines indicate the one-dimensional allowed regions for the two parameters at the same C.L. as the 2D plot.Figure 7
Allowed regions in thesin22θ13,Δmee2 space at a 68.27%, 95.45%, and 99.73% C.L. from our Daya Bay analysis. The best fit is marked by a “+” in the central plot, where the published Daya Bay contours [15] are also included. The marginalized Δχ2 over the oscillation parameters are also shown in the upper and right plots, including the rate-only result (dot-dashed line) for the mixing angle, and the 99.99% C.L. line (solid dark blue).Although the best fit is very well recovered, our contours are slightly wider than the published ones (dotted line in Figure7) towards the higher sin22θ13 values, and shorter towards the lower Δmee2. Despite our efforts to reproduce the full covariance matrix for this measurement, several manipulations had to be implemented in order to rebin the matrix and guarantee its positive definiteness, which may have introduced distortions. Nonetheless, we consider that our result captures the main features of the analysis and gives a good approximation to the confidence regions for the parameters. A cross-check calculation using a χ2 with pull terms produced contours with similar characteristics.


### 3.2.2. RENO

Following [25], we define the χ2 statistic for the RENO spectral analysis as
(12)χ2=∑i=1NbinsOiF/N−TiF/N2UiF/N+∑d=N,Fbdσbkdd2+eσscale2+∈σeff2.In this case,OiF/N is the ratio of the observed IBD candidate events at the far detector over those observed at the near detector for the ith energy bin; TiF/N=TiF/Nbd,e,∈,θ13,Δmee2 is the corresponding ratio of expected events, and UiF/N is the statistical uncertainty associated with OiF/N.The best fit found after minimizing theχ2 statistic with respect to the oscillation parameters, marginalizing over bd, e, and ∈ (the pull term parameters), is sin22θ13=0.083−0.012+0.010 and Δmee2=2.64−0.27+0.21×10−3eV2 at 1σ C.L., with χmin2/NDF=20.92/25. The results of this analysis are presented in Figure 8, where the allowed regions in the sin22θ13,Δmee2 parameter space are shown, together with the best-fit point. As in the case for the Daya Bay analysis, the 1-dimensional marginalized distributions are also shown for each of the oscillation parameters, and the result obtained from the rate-only analysis (dashed line in top panel) for sin22θ13 is included for comparison. We have also included the contours and best fit obtained by RENO [25] in Figure 8 (dotted line) which show good agreement with our analysis.Figure 8
Allowed regions in the (sin22θ13,Δmee2) space at a 68.27%, 95.45%, and 99.73% C.L. from our RENO analysis. The best fit is marked by a “+” in the central plot, where the published RENO contours [25] are also included. The marginalized Δχ2 over the oscillation parameters are also shown in the upper and right plots, including the rate-only result (dot-dashed line) for the mixing angle, and the 99.99% C.L. line (solid dark blue).The bottom-right plot in Figure3 compares the FD data to the no oscillations and best-fit predictions obtained from the measured spectrum at the ND. The agreement between the data and the best-fit prediction is very good and demonstrates that the near-to-far ratio technique used in the RENO analysis is insensitive to the presence of the 5 MeV bump.


## 3.2.1. Daya Bay

For Daya Bay, we follow Ref. [43] and define
(10)χ2=∑i,jNiobs−NiexpTVij−1Njobs−Njexp,where the indices i and j run over 156 bins corresponding to the concatenation of the prompt energy distributions of the 6 ADs, with 26 bins each. NiobsNiexp is the number of observed (expected) events in the ith energy bin, and Vij are the elements of the total covariance matrix expressed in the same binning. Niexp depends on the oscillation parameters sin22θ13 and Δmee2. The covariance matrix, V=Vstat+Vsyst, contains all sources of statistic and systematic errors affecting the experiment. The systematic error component was calculated from the full correlation matrix reported in [39], including the signal, background, and reactor core errors, and the total systematic uncertainties presented in Figure 2 of Ref. [40]. This assumption proved to be a reasonable approximation to the total systematic errors and correlations in the data set used for this analysis.The systematic error correlation matrix in [39] is a 222×222 matrix, since the energy spectrum in each AD used therein has 37 (nonuniform) bins in the energy range 0.7-12 MeV. However, the Daya Bay data set considered here used energy distributions with 26 bins, in the same energy range, hence having different bin boundaries (see Figure 5). In order to cast this correlation matrix in the form of a 156×156 matrix, we implemented a rebinning procedure based on the diagonal blocks (37×37 bins per AD) of the original matrix in which we sampled energy distributions consistent with the original matrix [44], and for each one, sampled 106 energy values which were then filled into a rebinned histogram (26×26 bins per AD). The resulting 1000 rebinned distributions were used to recalculate the correlation matrix with the desired binning. The full 6 AD correlation matrix, ρsyst, was constructed assuming a correlation among the 6 detectors encoded in a 6×6 matrix (see middle plot in Figure 6) which was adjusted so that the overall features of the original matrix could be reproduced. Our rebinned correlation matrix is shown in the top panel of Figure 6. The 156×156 elements of the total systematic covariance matrix Vsyst were then computed as
(11)Vijsyst=ρijsyst×Niexpσimod26×Njexpσjmod26(no summation over repeated indices) where the σk, k=1,⋯,26, are the total fractional systematic uncertainties, shown in the bottom panel of Figure 6, extrapolated from [40] up to 12 MeV. Despite being a rough approximation to the true error matrix used by the collaboration, as we will see, our results agree reasonably well with those reported by Daya Bay.Figure 5
Energy range and bin boundaries for the Daya Bay analysis in [39] (black lines and numbers), and for the analysis presented in this work (blue lines and numbers).Figure 6
(a) Full systematic error (signal, background, and reactor cores) correlation matrix used in the Daya Bay analysis with 26 bins per AD (signal, backgrounds, and reactor cores). Each square block represents the correlation among the bins of two ADs. (b) Correlation matrix for the 6 ADs. (c) Total fractional systematic uncertainties in the prompt energy distribution, assumed equal for all ADs.
(a)
(b)
(c)We found that the oscillation parameters which best reproduce the data aresin22θ13=0.091−0.009+0.012 and Δmee2=2.60−0.22+0.18×10−3eV2 at 1σ C.L., with χmin2/NDF=59.91/154. This result is shown in Figure 7, together with the 68.27%, 95.45%, and 99.73% C.L. allowed regions for the oscillation parameters space. In the upper (right) panel of Figure 7, we show the Δχ2 marginalization over Δmee2sin22θ13, where the minimum of the curve Δχ2=0 points to the best-fit value. We have included here the result of the rate-only analysis (dash-dotted line in the top panel) for comparison purposes. In the marginalization plots, the horizontal (vertical) lines indicate the one-dimensional allowed regions for the two parameters at the same C.L. as the 2D plot.Figure 7
Allowed regions in thesin22θ13,Δmee2 space at a 68.27%, 95.45%, and 99.73% C.L. from our Daya Bay analysis. The best fit is marked by a “+” in the central plot, where the published Daya Bay contours [15] are also included. The marginalized Δχ2 over the oscillation parameters are also shown in the upper and right plots, including the rate-only result (dot-dashed line) for the mixing angle, and the 99.99% C.L. line (solid dark blue).Although the best fit is very well recovered, our contours are slightly wider than the published ones (dotted line in Figure7) towards the higher sin22θ13 values, and shorter towards the lower Δmee2. Despite our efforts to reproduce the full covariance matrix for this measurement, several manipulations had to be implemented in order to rebin the matrix and guarantee its positive definiteness, which may have introduced distortions. Nonetheless, we consider that our result captures the main features of the analysis and gives a good approximation to the confidence regions for the parameters. A cross-check calculation using a χ2 with pull terms produced contours with similar characteristics.


## 3.2.2. RENO

Following [25], we define the χ2 statistic for the RENO spectral analysis as
(12)χ2=∑i=1NbinsOiF/N−TiF/N2UiF/N+∑d=N,Fbdσbkdd2+eσscale2+∈σeff2.In this case,OiF/N is the ratio of the observed IBD candidate events at the far detector over those observed at the near detector for the ith energy bin; TiF/N=TiF/Nbd,e,∈,θ13,Δmee2 is the corresponding ratio of expected events, and UiF/N is the statistical uncertainty associated with OiF/N.The best fit found after minimizing theχ2 statistic with respect to the oscillation parameters, marginalizing over bd, e, and ∈ (the pull term parameters), is sin22θ13=0.083−0.012+0.010 and Δmee2=2.64−0.27+0.21×10−3eV2 at 1σ C.L., with χmin2/NDF=20.92/25. The results of this analysis are presented in Figure 8, where the allowed regions in the sin22θ13,Δmee2 parameter space are shown, together with the best-fit point. As in the case for the Daya Bay analysis, the 1-dimensional marginalized distributions are also shown for each of the oscillation parameters, and the result obtained from the rate-only analysis (dashed line in top panel) for sin22θ13 is included for comparison. We have also included the contours and best fit obtained by RENO [25] in Figure 8 (dotted line) which show good agreement with our analysis.Figure 8
Allowed regions in the (sin22θ13,Δmee2) space at a 68.27%, 95.45%, and 99.73% C.L. from our RENO analysis. The best fit is marked by a “+” in the central plot, where the published RENO contours [25] are also included. The marginalized Δχ2 over the oscillation parameters are also shown in the upper and right plots, including the rate-only result (dot-dashed line) for the mixing angle, and the 99.99% C.L. line (solid dark blue).The bottom-right plot in Figure3 compares the FD data to the no oscillations and best-fit predictions obtained from the measured spectrum at the ND. The agreement between the data and the best-fit prediction is very good and demonstrates that the near-to-far ratio technique used in the RENO analysis is insensitive to the presence of the 5 MeV bump.


## 3.3. Combined Analysis

Finally, we performed a combined analysis of the two data sets by considering theχ2 statistic
(13)χcomb2=ΔχDayaBay2+ΔχRENO2.This definition will answer the specific question: how probable is it that both experimental results come from the same underlying oscillation model? [45]. Both the rate-only and the spectral analysis were performed using the corresponding Δχ2=χ2−χmin2 statistic appropriate for each case.For the combined rate-only analysis, we found that the data is best described with a value ofsin22θ13=0.088−0.006+0.008. This result is shown in Figure 9, where Δχ2=χcomb2−χcombmin2 is plotted as a function of sin22θ13 (solid purple line). We also plot here the rate-only results from the independent analyses of Daya Bay (dotted line) and RENO (dash-dotted line). Horizontal colored lines are drawn to mark the allowed intervals for the parameter at 1-4σ C.L., which are smaller for the combined analysis, indicating the expected enhancement in significance.Figure 9
ObtainedΔχ2 distribution for sin22θ13 from the Daya Bay+RENO rate-only combined analysis (full purple line), together with the results from Daya Bay (dotted line) and RENO (dash-dotted line), separately.For the spectral combined analysis, theχcomb2 is also built as the sum of the corresponding statistics used for Daya Bay and RENO, as in (13). The minimization of such a function gives the values of the oscillation parameters which produce the best fit to the data: sin22θ13=0.087−0.008+0.007, Δmee2=2.59−0.14+0.15×10−3eV21σ, with χcombmin2/NDF=0.18/2. Together with the best-fit point, Figure 10 shows the 68.27%, 95.45%, and 99.73% C.L. allowed regions in the studied parameter space.Figure 10
2D regions (1, 2, and 3σ) and 1D intervals (1, 2, 3, and 4σ) allowed for the oscillation parameters obtained from the combination of Daya Bay and RENO data, using the spectral analysis. The best fit is marked by a “+” in the central plot. The results from the independent analyses (RENO dot-dashed; Daya Bay dotted) are shown for comparison.We also show in Figure10 (top and right panels) the 1-dimensional Δχ2 distributions for each oscillation parameter, marginalizing over the other one, where we have included the results obtained separately for Daya Bay (dotted line) and RENO (dash-dotted line). Clearly, the combined result (solid line) produces smaller intervals for sin22θ13 and Δmee2.Finally, using the prescription described in Ref. [45], we evaluate the compatibility between the two data sets, by calculating the parameter goodness (PG) as the χ-squared probability Prob(χcombmin2; Pc), where Pc=2 is the number of parameters coupling the two data sets. In this case, we obtain a compatibility of PG=91.5%. We note that the small discrepancies between our results for Daya Bay and the published ones may be leading us to a different compatibility level from what could be obtained with the official analysis results. Our wider contour in sin22θ13 tends to increase the compatibility, while being smaller along Δmee2 tends to reduce it. However, with regard to the question at the beginning of this section, the compatibility found here allows us to state that both the Daya Bay and RENO data considered in this work are well described by the same neutrino oscillation model, represented by (7), with the oscillation parameters found in the combined analysis.


## 4. Conclusions

We have studied the neutrino oscillation analyses from two particular data taking periods of the Daya Bay and RENO experiments. We reproduced reasonably well the published rate-only and spectral analysis results from both experiments, obtaining, for the spectral analysis:(Daya Bay)sin22θ13=0.091−0.012+0.009,Δmee2=2.60−0.22+0.18×10−3eV2,(RENO)sin22θ13=0.083−0.012+0.010,Δmee2=2.64−0.27+0.21×10−3eV2.The spectral analysis of Daya Bay was the more challenging, considering our choice to use the full systematic error covariance matrix in 26 prompt positron energy bins in the definition of theχ2 statistic. This required the implementation of a statistical method to rebin the correlation matrix found in a more recent publication by the collaboration. As a cross-check, we obtained very similar contours from a definition of the Daya Bay χ2 statistic using pull terms. We were able to reproduce very closely all the results of the RENO spectral analysis and verify that the near/far ratio technique makes the results insensitive to the presence of the 5 MeV bump.A combined spectral analysis was carried out by defining aχ2 statistic as the sum of the parameter goodness (PG) Δχ2 for each data set, and extracting confidence regions around its minimum. We found that the values that best fit the data are
(Daya Bay+RENO)sin22θ13=0.087−0.008+0.007,Δmee2=2.59−0.14+0.15×10−3eV2,at 1σ C.L. The combined analysis provided more restricted allowed regions for the oscillation parameters, compared against the results from the two experiments separately, as expected, with an increase in the precision of the oscillation parameters from 30 to 40%. Furthermore, we found the data sets considered here to be compatible at the 91.5% level according to our analyses, despite small discrepancies with our result and the one published by Daya Bay.



---

*Source: 8526034-2020-01-22.xml*