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 parameters and the effective mass splitting 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 and 40% with respect to the individual analyses considered.

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, , weakly interacting particles which are found to exist in three different flavors: electron neutrinos , 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 as where are the elements of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix, which depend on the mixing angles and a CP-violating phase [4]. The PMNS matrix may also depend on two additional Majorana phases , 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 in vacuum is given by [5, 6] where is the neutrino energy and are the differences of the squared masses of the definite-mass states and .

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 , by observing the disappearance of reactor antineutrinos 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 , the neutrino mass ordering, and the octant of the angle (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 , 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 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 with energies below 1.3 MeV [26].

Reactor antineutrinos with energies 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 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 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 [2731] 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].

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 [34], and uses the 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 , 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 ’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.

In the target volume, a large amount of freely moving protons may interact with the antineutrinos coming from the reactors, producing positrons , which then annihilate with surrounding electrons generating two gamma rays (prompt signal). In addition, in the IBD process, a neutron 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 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, , 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.

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.

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 reactor and arrives at a specific detector . 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 is the thermal power of reactor , is the mass of the fiducial volume of detector , and is the baseline distance between reactor and detector ; is a normalization constant making the sum .

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.

The neutrino energy for each simulated event is then assigned using the well-known relation [13] where is the average energy taken by the neutron (∼10 keV). A small multiplicative correction factor was applied to 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 , 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 Section 2.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 , prompt positron energy , and baseline are assigned as follows: first, a baseline 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 () is added to approximately incorporate the reactor and detector sizes. The prompt positron energy is then sampled from the distribution corresponding to the chosen detector (and corrected by multiplying it by ). This value is then used to calculate the true neutrino energy 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 probability in Equation (2) is well approximated by with , and is given in eV2, (in m) is the distance between the reactor and the detector, and (in MeV) is the neutrino energy. Given that , the oscillation is mainly driven by , and Equation (6) naturally leads to the definition of an effective squared-mass difference such that [24], so that the survival probability in Equation (6) can be reduced to

For the purpose of reproducing the Daya Bay results, in what follows, we have used , , for the spectral analysis, and also , for the rate-only analysis, as in Ref. [24]. For the reproduction of the RENO results, we used , and , for the spectral analysis, and in addition for the rate-only analysis, as in Ref. [25].

It has been pointed out that the definition of used here, besides being dependent, is discontinuous at 0.5 km/MeV [42], and better definition can be considered, such as the weighted average of and . 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 Tables 1 and 2, and the measured values of the oscillation parameters , , and , we extract the value of . Since the shape of the spectrum is not considered, no information about is obtained from this analysis.

For Daya Bay, we follow Ref. [15] and define our statistic as

Here, is the number of IBD-observed events in the th detector after background subtraction, is the background rate, and is the predicted number of events considering neutrino oscillations through Equation (7), in the total DAQ live time (Table 1). The quantity is the fraction of events produced in reactor which contribute to detector , considering the travel distance and the neutrino flux, which corresponds precisely with the baseline probability of Figure 4. The in (8) is penalized by the inclusion of 18 pull terms , , and (with ), characterizing the systematic errors affecting the measurement. As reported by the collaboration in Ref. [15], (0.8%) and (0.2%) are the uncorrelated reactor and detector uncertainties, respectively, and 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 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 statistic. Minimization of the marginalized gives the best-fit value 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 statistic as where is the number of observed events after background subtraction for the near () and far () detectors in the total DAQ live time (Table 2); is the number of expected events in detector coming from reactor , including the detection efficiency and the effect of oscillations. Here, is a freely varying normalization factor, , from Table 2, are the background uncertainties associated with the pull term parameters , and (0.9%) and (0.2%) are the uncorrelated reactor and detector uncertainties, associated with the pull term parameters and , respectively. Using a similar minimization procedure to the one used for the Daya Bay analysis, marginalizing over and the pull term parameters, we obtain 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 along with . The best fit is found by minimizing a suitable statistic defined over a uniformly spaced grid in the vs. space.

3.2.1. Daya Bay

For Daya Bay, we follow Ref. [43] and define where the indices and run over 156 bins corresponding to the concatenation of the prompt energy distributions of the 6 ADs, with 26 bins each. is the number of observed (expected) events in the th energy bin, and are the elements of the total covariance matrix expressed in the same binning. depends on the oscillation parameters and . The covariance matrix, , 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 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 matrix, we implemented a rebinning procedure based on the diagonal blocks ( 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 ( 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, , was constructed assuming a correlation among the 6 detectors encoded in a 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 elements of the total systematic covariance matrix were then computed as (no summation over repeated indices) where the , , 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.

We found that the oscillation parameters which best reproduce the data are and at 1σ C.L., with . 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 marginalization over , where the minimum of the curve 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 (dotted line in Figure 7) towards the higher values, and shorter towards the lower . 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 with pull terms produced contours with similar characteristics.

3.2.2. RENO

Following [25], we define the statistic for the RENO spectral analysis as

In this case, is the ratio of the observed IBD candidate events at the far detector over those observed at the near detector for the th energy bin; is the corresponding ratio of expected events, and is the statistical uncertainty associated with .

The best fit found after minimizing the statistic with respect to the oscillation parameters, marginalizing over , , and (the pull term parameters), is and at 1σ C.L., with . The results of this analysis are presented in Figure 8, where the allowed regions in the 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 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.

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.

3.3. Combined Analysis

Finally, we performed a combined analysis of the two data sets by considering the 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 statistic appropriate for each case.

For the combined rate-only analysis, we found that the data is best described with a value of . This result is shown in Figure 9, where is plotted as a function of (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.

For the spectral combined analysis, the 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: , , with . 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 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 and .

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(; ), where is the number of parameters coupling the two data sets. In this case, we obtain a compatibility of . 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 tends to increase the compatibility, while being smaller along 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:

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 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 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 statistic as the sum of the parameter goodness (PG) for each data set, and extracting confidence regions around its minimum. We found that the values that best fit the data are 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.

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].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

M.A.A. and D.J.P.-T. thank the Organizing Committee of the VIII Encuentro Regional de Ciencias Físicas 2018, Barranquilla, Colombia, in which preliminary results of this work where presented. M.A.A. also thanks the Instituto de Ciencias Nucleares for their hospitality during the partial realization of this work. A.A.A.A. acknowledges the hospitality of the Universidad del Atlántico during the time he spent working at the Physics Department. This work was supported by the Universidad del Atlántico through the grant “Convocatoria Interna Impacto Caribe” No. CB71-CIC2014 and by the Consejo Nacional de Ciencia y Tecnología (CONACyT, México) through SNI (Sistema Nacional de Investigadores).