Neutrino oscillation analysis of 217 live-days of Daya Bay and 500 live-days of RENO

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 parameters sin 2θ13 and the effective mass splitting ∆m 2 ee 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 χ statistic for each case. Finally, we performed a combined spectral analysis and extract tighter constraints on the parameters, with an improved precision between 30-40% with respect of the individual analyses considered.


I. 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 non-zero 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 ν k (k = 1, 2, 3) as 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] P να→ν β = δ αβ − 4 3 k>j Re U * αk U βk U αj U * βj sin 2 ∆m 2 where E is the neutrino energy and ∆m 2 kj ≡ m 2 k −m 2 j 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 * marioacero@mail.uniatlantico.edu.co † alexis@nucleares.unam.mx ‡ djosepolo@mail.uniatlantico.edu.co [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 non-zero value, and opening the door to studying CP violation in the neutrino sector. In later years, neutrino oscillations 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] and [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.

II. ANTINEUTRINOS FROM NUCLEAR REACTORS
Typical commercial pressurized water reactors (PWR) are copious sources ofν e s, originating primarily in the beta decay of unstable fission products of the fissile isotopes 235 U, 238 U, 239 Pu, and 241 Pu, as well as from neutron capture in 238 U. On average, each fission releases roughly 200 MeV of energy and produces 6 antineutri-nos with energies below 10 MeV. Since the typical decay chain of the fission products has three consecutive beta decays, about 2 × 10 20ν e per second per Giga-Watt of thermal power (GW th ) [6] are isotropically emitted from the reactor core. The neutron capture contribution occurs at a smaller rate (0.6 per fission) and producesν e s with energies below 1.3 MeV [26].
Reactor antineutrinos with energies > 1.8 MeV can be detected via the inverse beta decay (IBD) process 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ν e s 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 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 non-standard particle physics (see [27][28][29][30][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.

A. 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 GW th 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 (AD). In the data taking period used in this work, the detectors were distributed in three experimental halls (EH1, EH2, 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 [33].
The 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 [37], and uses theν e from six PWRs arranged in a line along the coast. The reactors produce a total of 16.4 GW th . RENO uses two detectors (Near Detector -ND-and Far Detector -FD-) to observe the producedν e s, as displayed in Figure 1(b), where the red circles represent the six reactors, and the blue circles represents 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, Eq. (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 [33] for DB and [36] for RENO.
In the target volume, a large amount of freely moving protons (p) 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  [34]) and RENO (bottom [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.
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]).

B. Input to our studies
The prompt reconstructed energy, E p , distributions used in the analyses are shown in Figures 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 nooscillation 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 II C). 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.
In order to normalize the event rates and energy distributions, we collected the information from tables found in Refs. [24] and [25], which we have summarized here in Table I and Table II, for Daya Bay and RENO, respectively. In these tables we have added 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. [42], the total systematic errors from figure 2 in Ref. [44], and used this information to construct the full covariance matrix, as will be described in section III below.

C. 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 reactor i and arrives at a specific detector j. 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: where P th r is the thermal power of reactor r, M d is the mass of the fiducial volume of detector d, and L rd is the baseline distance between reactor r and detector d; N is a normalization constant making the sum rd w d r = 1. The baseline lengths and detector fiducial volume masses were extracted from Ref. [33] for the Daya Bay analysis, and from Ref. [37] 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 reaching a particular detector, was produced in a specific reactor. As expected, near detectors have larger probabilities than far detectors, since the closer the detector is to a reactor, the more antineutrinos are detected.
The neutrino energy for each simulated event is then assigned using the well-known relation [13] where E n is the average energy taken by the neutron (∼ 10 keV). A small multiplicative correction factor f = 1.002 was applied to E p (E p → E p × f ) in our simulation to better reproduce the results. 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 section II B, Refs. [24] and [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 Eq. (7) to a sample of 10 7 simulated neutrino IBD events whose true energy E ν , prompt positron energy E p , 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 (σ = 1 m) is added to approximately incorporate the reactor and detector sizes. The prompt positron energy is then sampled from the E p distribution corresponding to the chosen detector (and corrected by multiplying it by f ). This E p value is then used to calculate the true neutrino energy E ν using Eq. (5).
It has been pointed out that the definition of ∆m 2 ee used here, besides being L/E dependent, is discontinuous at 0.5 km/MeV [39], and better definition can be considered, such as the weighed average of ∆m 2 31 and ∆m 2 32 . 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 not 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.

A. Rate-only analysis
Using only the information of the total event rates reported in Table I and Table II, and the measured values of the oscillation parameters sin 2 2θ 12 , ∆m 2 21 and ∆m 2 32 , we extract the value of sin 2 2θ 13 . Since the shape of the spectrum is not considered, no information about ∆m 2 ee is obtained from this analysis.
For Daya Bay, we follow Ref. [15] and define our χ 2 statistic as: Here, M d is the number of IBD observed events in the dth detector after background subtraction, B d is the background rate, and T d is the predicted number of events considering neutrino oscillations through Eq. (7), in the total DAQ live time (Table I). The quantity ω d r 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 interval 0 < sin 2 2θ 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 sin 2 2θ 13 = 0.090 +0.010 −0.009 at 1σ C.L.
A similar rate-only analysis was performed to the RENO data, considering the three-neutrino oscillation model described by Eq. (7). In this case, we follow Ref. [16] and define the χ 2 statistic as where N d obs is the number of observed events after background subtraction for the near (N ) and far (F ) detectors in the total DAQ live time (Table II); N d,r exp 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, σ d b , from Table II are the background uncertainties associated to the pull term parameters b d , and σ r (0.9%) and σ ξ d (0.2%) are the uncorrelated reactor and detector uncertainties, associated to the pull term parameters f r 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 sin 2 2θ 13 = 0.088 +0.010 −0.013 at 1σ as the value which best fit to the RENO data.

B. 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 ∆m 2 ee along with sin 2 2θ 13 . The best fit is found by minimizing a suitable χ 2 statistic defined over a uniformly spaced 100×100 grid in the sin 2 2θ 13 vs ∆m 2 ee space.

Daya Bay
For Daya Bay we follow Ref. [41] and define ) is the number of observed (expected) events in the i-th energy bin, and V ij are the elements of the total covariance matrix expressed in the same binning. N exp i depends on the oscillation parameters sin 2 2θ 13 and ∆m 2 ee . The covariance matrix, V = V stat + V syst , contains all sources of statistic and systematic errors affecting the experiment. The systematic error component was calculated from the full correlation matrix reported in [42], including the signal, background and reactor core errors, and the total systematic uncertainties presented in figure 2 of Ref. [44]. 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 [42] is a 222 × 222 matrix, since the energy spectrum in each AD used therein has 37 (non-uniform) 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 Fig. 5). In order to cast this correlation matrix in the form 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 [43], and for each one, sampled 10 6 energy values which were then filled into a re-binned histogram (26 × 26 bins per AD). The resulting 1000 re-binned distributions were used to re-calculate 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 top-right plot in Figure 6) which was adjusted so that the overall features of the original matrix could be reproduced. Our re-binned correlation matrix is shown in the left panel of Figure 6. The 156×156 elements of the total systematic covariance matrix V syst were then computed as (no summation over repeated indices) where the σ k , k = 1, . . . , 26, are the total fractional systematic uncertainties, shown in the bottom-right panel of Figure 6, extrapolated from [44] 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  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 ∆m 2 ee (sin 2 2θ 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.
Although the best fit is very well recovered, our contours are slightly wider than the published ones towards (dotted line in Fig. 7) the higher sin 2 2θ 13 values, and shorter towards the lower ∆m 2 ee . Despite our efforts to reproduce the full covariance matrix for this measurement, several manipulations had to be implemented in order to re-bin 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.

RENO
Following [25] we define the χ 2 statistic for the RENO spectral analysis as In this case, O  Figure 8, where the allowed regions in the (sin 2 2θ 13 , ∆m 2 ee ) 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 sin 2 2θ 13 is included for comparison. We have also included the contours and best fit obtained by RENO [25] in Fig. 8 (dotted line) which show good agreement with our analysis.
The bottom-right plot in Figure 3 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.

C. Combined Analysis
Finally, we performed a combined analysis of the two data sets by considering the χ 2 statistic 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 −χ 2 min statistic appropriate for each case.
For the combined rate-only analysis we found that the data is best described with a value of sin 2 2θ 13 = 0.088 +0.008 −0.006 . This result is shown in Figure 9, where ∆χ 2 = χ 2 comb − χ 2 comb min is plotted as a function of sin 2 2θ 13 (solid purple line). We also plot here the rateonly 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.
For the spectral combined analysis, the χ 2 comb 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: sin 2 2θ 13 = 0.089 +0.007 −0.008 , ∆m 2 ee = 2.59 +0.15 −0.14 × 10 −3 eV 2 (1σ), with χ 2 comb min /NDF = 0.31/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.
We also show in Figure 10 (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 sin 2 2θ 13 and ∆m 2 ee . 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(χ 2 comb min ; P c ), where P c = 2, is the number of parameters coupling the two data sets. In this case we obtain a compatibility of PG = 86%. We note that the small discrepancies between our results for Daya Bay and the published ones, are pulling the compatibility to a lower value than could otherwise be expected. However, with regards 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.

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 analyses results from both experiments, obtaining, for the spectral analysis:  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 re-bin 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 sin 2 2θ 13 = 0.089 +0.007 −0.008 , ∆m 2 ee = 2.59 +0.15 −0.14 × 10 −3 eV 2 [Daya Bay + RENO] 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-40%. Furthermore, we found the data sets considered here to be compatible at the 86% level according to our analyses, despite small discrepancies with our result and the one published by Daya Bay.

DATA AVAILABILITY
The data used to support the findings of this analyses are included within the article for better readability, and are also openly accessible in [24,25].