Search for a time-varying electron antineutrino signal at Daya Bay

A search for a time-varying $\bar{\nu}_{e}$ signal was performed with 621 days of data acquired by the Daya Bay Reactor Neutrino Experiment over 704 calendar days. The time spectrum of the measured $\overline{\nu}_e$ flux normalized to its prediction was analyzed with a Lomb-Scargle periodogram, which yielded no significant signal for periods ranging from 2 hours to nearly 2 years. The normalized time spectrum was also fit for a sidereal modulation under the Standard Model extension (SME) framework to search for Lorentz and CPT violation (LV-CPTV). Limits were obtained for all six flavor pairs $\bar{e}\bar{\mu}$, $\bar{e}\bar{\tau}$, $\bar{\mu}\bar{\tau}$, $\bar{e}\bar{e},\bar{\mu}\bar{\mu}$ and $\bar{\tau}\bar{\tau}$ by fitting them one at a time, constituting the first experimental constraints on the latter three. Daya Bay's high statistics and unique layout of multiple directions from three pairs of reactors to three experimental halls allowed the simultaneous constraint of individual SME LV-CPTV coefficients without assuming others contribute negligibly, a first for a neutrino experiment.

A search for a time-varyingνe signal was performed with 621 days of data acquired by the Daya Bay Reactor Neutrino Experiment over 704 calendar days. The time spectrum of the measured νe flux normalized to its prediction was analyzed with a Lomb-Scargle periodogram, which yielded no significant signal for periods ranging from two hours to nearly two years. The normalized time spectrum was also fit for a sidereal modulation under the Standard-Model Extension (SME) framework to search for Lorentz and CPT violation (LV-CPTV). Limits were obtained for all six flavor pairs eμ,ēτ ,μτ ,ēē,μμ andττ by fitting them one at a time, constituting the first experimental constraints on the latter three. Daya Bay's high statistics and unique layout of multiple directions from three pairs of reactors to three experimental halls allowed the simultaneous constraint of individual SME LV-CPTV coefficients without assuming others contribute negligibly, a first for a neutrino experiment.

I. INTRODUCTION
Some scenarios of physics beyond the Standard Model (SM) predict a time-varying probability of neutrino oscillation. Among these are models in which ultra-light scalar dark matter couples to neutrinos, inducing periodic variations in the mass splittings and mixing angles [1,2]. Other models involve Lorentz symmetry violation (LV), which is suggested as a signature of Planck scale phenomenology [3][4][5][6] and which could be accompanied by CPT violation (CPTV) [4,7].
The Standard-Model Extension (SME) [8][9][10] was introduced as an effective theory that maintains the usual gauge structure and properties of the SM such as renormalizability, but adds all the possible terms constructed with SM fields that introduce Lorentz symmetry breaking. By predicting a set of testable signatures in various areas of physics, it provides a connection between experimental research and more fundamental theories extending to the Planck scale.
In the neutrino sector, the violation of rotationsymmetry in the SME causes deviations from standard oscillation probabilities derived from the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [11] that depend on propagation direction. This would produce a timevarying neutrino oscillation probability associated with the Earth's orbital and rotational movement relative to the fixed stars, and therefore a period of a sidereal day (23 h 56 min 4.09 sec). Accordingly, a sidereal time dependence has been sought in the oscillation probability of * Now at Department of Chemistry and Chemical Technology, Bronx Community College, Bronx, New York 10453 accelerator neutrinos [12][13][14][15][16][17], atmospheric neutrinos [18] and reactor neutrinos [19]. The SME also predicts deviations from the standard L/E oscillation behavior. The oscillated neutrino energy spectrum of atmospheric neutrinos has been examined for such a distortion both in the Super Kamiokande [20] and IceCube [21] experiments. No positive LV or CPTV signal has yet been observed, and neutrino oscillation experiments have set some of the most stringent limits on the violation of these fundamental symmetries of nature, down to the level of 10 −28 [21].
The Daya Bay reactor neutrino experiment [22] has recently produced the most precise measurements of reactor electron antineutrino (ν e ) disappearance at short baselines [23][24][25][26][27][28]. In neutrino oscillation experiments, time-dependent LV-CPTV effects are amplified with distance, and Daya Bay's baselines (< 2 km) are relatively short compared to many other neutrino oscillation experiments. However, Daya Bay has accumulated the largest sample of reactor ν e 's to date. Moreover, it has a unique experimental layout comprising different well-known neutrino propagation directions. Both of these factors make it an excellent experiment to search for a time-varyingν e signal and LV-CPTV effects.
This letter first describes a generic search for an unpredicted periodicity in theν e rates measured at Daya Bay using the Lomb-Scargle method [29]. This analysis, which yields no positive results, has the potential to identify the presence of an unexpected time-variant source of ν e 's. The letter then presents a targeted search for a sidereal time modulation in the context of the SME, producing limits on the individual coefficients that characterize the theory. The Daya Bay reactor complex consists of three nuclear power plants (Daya Bay, Ling Ao, Ling Ao II), each with two reactors. The emitted ν e flux is sampled in eight identically designed antineutrino detectors (ADs) located in three experimental halls (EHs), as shown in Fig. 1. Each AD is filled with 20 tons of gadolinium-doped liquid scintillator enclosed by 22 tons of undoped liquid scintillator and 40 tons of mineral oil. Scintillation light is detected by 192 photomultiplier tubes (PMTs). The ADs are immersed in pure water pools that are instrumented with PMTs, providing shielding and serving as cosmogenic muon detectors. Candidate events independently trigger each detector and are read out by custom front-end electronics. For the ADs, a readout window of 1.2 µs of data is initiated when the number of PMTs with an above-threshold signal is greater than 45 or the synchronous analog sum of charge output by the PMTs is larger than a value corresponding to ∼0.4 MeV [22]. The clock for the readout electronics and trigger systems runs at 40 MHz, and is synchronized to a global 10 MHz signal generated by a rubidium oscillator further synchronized to absolute Coordinated Universal Time (UTC) with a Global Positioning System (GPS) receiver. Further information about the Daya Bay experiment can be found in Ref. [22].

B. Antineutrino Signal and Backgrounds
The data set used in this study corresponds to a total exposure of 621 days distributed over 704 solar days (705.5 sidereal days), from December 24, 2011 to November 27, 2013. Data-taking began with 6 ADs and continued for 217 days, pausing from July 28, 2012 to October 19, 2012 for the installation of the final two ADs, one in EH2 and the other in EH3. Each physics run lasted as long as 72 hours. A three-hour interruption of normal data acquisition occurred almost every Friday for calibrating the detectors.
Electron antineutrinos were detected via the inverse beta decay (IBD) reaction, ν e + p → e + + n, where the energy loss of the e + in the scintillator and its subsequent annihilation provided a prompt scintillation light signal followed by a delayed light signal from the neutron capture on gadolinium. IBD candidates were selected by requiring prompt-delayed pairs to have specific energies (0.7 < E prompt < 12.0 MeV, 6.0 < E delayed < 12.0 MeV) and time separation (1 < ∆t < 200 µs). Selected events were also required not to have been preceded by a muon candidate. A multiplicity cut was applied to ensure that only isolated prompt-delayed pairs were selected. Two slightly different IBD selections based on these criteria were used in two independent analyses, which also estimated backgrounds differently. Distinct muon veto and multiplicity cut efficiencies were accurately assessed from muon and random background rates as a function of time, and applied in the estimation of the IBD rates. The timedependence of these efficiencies was negligible.
The total background amounted to less than 3% of the total IBD candidate samples and was dominated by accidental coincidences. Both analyses precisely determined this background hourly from the measured rates of uncorrelated signals and subtracted it from the IBD samples. One analysis also considered the small variations in the fast neutron and 9 Li/ 8 He correlated backgrounds, which were determined for the full period and then estimated hourly by scaling them with the measured muon rate. The same was done for the background caused by neutrons from the 241 Am-13 C calibration sources, but scaling with the hourly single neutron rate. The slight timedependence in these backgrounds was found to contribute < 0.01% of the total variation in IBD rates, and thus had a negligible impact on the results presented in this letter.
The ν e oscillation probability was determined as the ratio of the measured IBD rate to the predicted IBD rate assuming no oscillation. The measured IBD rate in each hall was determined hourly by dividing the number of background-subtracted IBD events by the data acquisition livetime, correcting for the loss of time caused by the muon veto and multiplicity cuts. This rate was then divided by the hourly expectation determined as in Ref. [31] but with a livetime-weighted linear interpola- tion of daily thermal power data, yielding the measured survival probability P (t). The time average of P (t) was normalized to the Lorentz-invariant three-neutrino survival probability P (0) measured by Daya Bay in Ref. [26] with the same data set, and P (0) was subtracted to give the residual survival probability R(t) ≡ P (t)−P (0) . This quantity is shown in Fig. 2 and Fig. 3, the former vs. real time in hourly bins, and the latter accumulated in 24 sidereal hour bins [30]. The origin for the sidereal time was set to local midnight on the 2018 Vernal equinox, a convention typically used by experiments searching for LV and CPTV in the sun-centered frame. An integer number of sidereal days was subtracted to obtain the closest foregoing time to the start of data-taking, yield-ing 2011/12/23 22:13:13.80705 UTC. This choice had no impact on the limits reported in Sections III and IV.

C. Uncertainties
Given the nature of this search, only uncertainties of quantities that varied over time were taken into account. This included statistical, reactor-related, and event selection uncertainties. The statistical uncertainty dominated the uncertainty in each EH, contributing at the level of 0.63%, 0.71%, and 1.26% of P (t) in each of the 24 time bins of Fig. 3 for EH1, EH2 and EH3, respectively.
The efficiency uncertainty was dominantly due to the delayed energy (E delayed ) cut, and was inferred using the estimated stability of the energy scale. The energy scale was calibrated during data collection using spallation neutrons, and was found to vary within 0.2% in all ADs [22]. Variations in the number of target protons, amount of neutrons produced by IBD interactions outside the target that diffused into the target, and neutron capture time were estimated with the relationship between the density ρ and temperature T of the gadolinium-doped liquid scintillator: ∆ρ = −9.05 × 10 −4 ∆T [32]. The expected 0.045% change in ρ based on the observed 0.5 • C variation in temperature was propagated to the uncertainties of these parameters. Uncertainties of all other selection efficiencies were less significant and conservatively inherited from the oscillation analysis [26]. The overall uncertainty of the selection efficiency in the three halls was estimated to be 0.09% of the survival probability P (t) for each bin of Fig. 3. When combining the data of individual ADs in the same hall, correlations were considered.
All the uncorrelated reactor-related uncertainties involved in the flux prediction, which included power, energy/fission, fission fraction, spent fuel and nonequilibrium corrections, totaled to 0.9% and were conservatively treated as time-dependent on a daily basis. The relative size of the reactor systematic with respect to the survival probability P (t) for each bin in Fig. 3 was 0.10%, 0.09%, and 0.08% for EH1, EH2, and EH3, respectively. Correlations between the predicted fluxes at the three halls had negligible impact to the analyses presented.
As discussed previously, background variation with time was found to be negligible.

III. ANALYSIS ON PERIODIC AMPLITUDES
A general search for a periodic signal within the measured residual survival probability was performed for each of the three experimental halls using the Lomb-Scargle (LS) periodogram [29], which is a widely used technique for detecting periodic signals in unevenlysampled data. A periodogram was derived for each panel in Fig. 2, spanning a frequency range from 5.9×10 −5 sidereal hour −1 to 0.5 sidereal hour −1 . The normalized LS power for a frequency f derived from N data points X j at specific times t j can be estimated as [29] with X ≡ N j=1 X j /N and τ defined by tan(4πf τ ) = N j=1 sin(4πf t j )/ N j=1 cos(4πf t j ). The normalization is accomplished by dividing by the total variance, σ 2 ≡ N j=1 (X j − X) 2 /(N − 1). The obtained values of σ 2 were 0.023, 0.032 and 0.112 for EH1, EH2 and EH3, respectively. The bottom panels of Fig. 4 show the resulting LS powers for each frequency in each hall.
It was noted in Ref. [33] that if the signal X j is purely white noise, then L(f ) follows an exponential probability distribution when normalized with σ 2 . Accordingly, the significance of a given LS power can be determined with a confidence level (CL) defined as (1 − e −L(f ) ) M , where M is the number of independent frequencies that are scanned. M is nearly equal to the number of data points M ≈ N in the case of even sampling, but is a priori unknown for unevenly distributed samples. To estimate this number, 10,000 Monte Carlo data sets with statistical fluctuations were analyzed. The highest LS power z in each data set was selected to construct a probability density function for each hall, which was then fit as [33,34]. The extracted values of M were 16588, 16245 and 16697 for EH1, EH2, and EH3, respectively, while N = 16913. The variations were caused by the different statistics of each hall as modeled in the simulated data sets, and had little impact on the CLs.
The resulting CL values for each frequency can be seen on the top panels of Fig. 4. Table I gives information about the highest LS power in each EH. It is noteworthy that none of the highest powers are common among the three halls. No significant evidence for a periodic signal was found.
The periodicity search was also performed with the Discrete Fourier Transform (DFT). Since this method does not account for uneven sampling, the gaps in data acquisition were handled by exploiting the linearity of the transform. The DFT was applied to the data, with the residual survival probability set to zero in all the gap bins (see Fig. 2). Many simulated data sets with the same gaps and no time-varying signal were also transformed and averaged for each hall, and then the results subtracted from the data. The impact of the subtraction was very small, which is expected due to the small number of missed hourly samples (typically, a few per week) relative to the total number of samples (ideally, 168 per week). The resulting power spectra were consistent with those obtained from the LS method.

IV. ANALYSIS ON LV-CPTV COEFFICIENTS
The data were also probed for an LV-CPTV signal under the SME. In this framework, the survival probability can be expressed as P νe→νe = P (0) + P (1) + P (2) + .... The first term P (0) = |S (0) eē | 2 is the mass-driven survival probability for ν e 's in the Lorentz-invariant case. P (1) is calculated as [35] where L is the baseline, (M (1) eē )cd are the so-called experimental factors, T ⊕ represents sidereal time and ω ⊕ = 2π/(1 sidereal day). The subscriptcd runs over theēē, µμ,ττ ,ēμ,ēτ andμτ flavor pairs. (C)cd, (A s )cd, (A c )cd, (B s )cd and (B c )cd are commonly referred to as the sidereal amplitudes, which are functions of a total of fourteen SME coefficients for each flavor pair, as well as neutrino energy and propagation direction. The complex relationship between the sidereal amplitudes and the individual coefficients, as well as other details concerning the SME, can be found in the appendix.
The goal of this analysis is to constrain the individual SME coefficients contained in the sidereal amplitudes. For Daya Bay's baselines and energies, P (2) is smaller than P (1) by a few orders of magnitude, and was consequently ignored together with higher order terms. The subtraction of P (0) in R(t) made the fit insensitive to the isotropic amplitude (C)cd, whose coefficients can be extracted by analyzing the time-independent energy and baseline dependencies of the oscillation probability. These effects have been constrained by atmospheric neutrino data [20,21] well beyond the reach of Daya Bay. Without (C)cd, a total of nine different coefficients are contained in the amplitudes (A s )cd, (A c )cd, (B s )cd and (B c )cd, as shown in Eq. (A8). The sum in P (1) over the six flavor pairs makes it unfeasible for a single experiment to simultaneously constrain the 8 × 6 = 48 parameters with one fit, given their degeneracies. Interplay between the terms could be disentangled by comparing results from experiments with different neutrino energies, directions, and flavors. Without a positive signal however, it is impossible to determine whether there are any correlations or cancellations among the terms within the sum. Accordingly, the standard practice of fitting each flavor pair at a time by setting the coefficients of the other pairs to zero was employed.
Up to now neutrino experiments [12][13][14][15][16][17][18][19] have reported limits on the four sidereal amplitudes (A s )cd, (A c )cd, (B s )cd and (B c )cd. Even when considering individual flavor pairs, the number of directiondependent parameters precluded these experiments from setting limits on individual coefficients, except through the method of fitting one coefficient at a time while arbitrarily setting all others to zero. With a unique configuration of multiple directions from three experimental sites to three pairs of nuclear reactors (see Fig. 1) and the separation into five energy bins described below, Daya Bay was able to completely disentangle the energy and direction dependencies in Eq. (A8) and to simultaneously constrain eight LV-CPTV coefficients: , chosen so as to contain a similar amount of statistics in each. This resulted in 15 independent data sets (3 EHs × 5 energy bins), whose residual survival probabilities R j (t) = P j (t) − P (0) j are shown in Fig. 5. Given that each EH sees ν e 's from the six reactor cores, these data sets were simultaneously fit with where f ij is the expected fraction of events from the ith reactor core in data set j, and P (1) ij is the oscillation probability of Eq. (2) for that particular ij combination.
The event fraction f ij was calculated as where F ij is the i-th core's time-integrated flux seen in the hall corresponding to data set j, determined from the reactor power and fission fraction information provided by the power plant [31] and including oscillation and inverse-square law effects. Accordingly, the χ 2 used in the fit is expressed as Here R = P (t) − P (0) is the measured residual survival probability of each data point, R fit is the SME prediction given by Eq. (3), and σ 2 R is the total error as described in Section II C. The energy spread in each of the five bins was taken into account when calculating P (1) ij and found to be unimportant. Prompt energy was converted to ν e energy using a response matrix [28]. The fit was performed assuming the normal neutrino mass ordering, a zero value for the CP -violating phase, and the values of the oscillation parameters reported in Ref. [36]. The first two choices, as well as the uncertainties of the oscillation parameters, were found to have a negligible impact on the results.
The two analyses obtained very similar best-fit values and limits, which are shown in Table II  As a probe of new physics, a model-independent search for a time-variation of the reactor ν e survival probability was performed with 621 days of Daya Bay data over a period of 704 calendar days. The Lomb-Scargle method yielded no significant evidence for a periodicity in the frequency range of 5.9×10 −5 sidereal hour −1 to 0.5 sidereal hour −1 . The survival probability measured at Daya Bay was also examined for a sidereal time dependence within the SME framework. Daya Bay's high statistics and multiple-baseline configuration allowed a complete disentangling of the energy and direction dependencies within the sidereal amplitudes, yielding the first simultaneous constraints of individual Lorentz-violating coefficients for a neutrino experiment. Limits were provided for theēμ,ēτ ,μτ ,ēē,μμ andττ flavor pairs, yielding the first experimental constraints for the latter three.
This appendix summarizes the details involved in calculating the SME prediction for the analysis presented in Section IV and lays out the relationship between the sidereal amplitudes and the individual coefficients. Supplementary materials providing necessary values for the reader to reproduce the results presented in this letter are included online [30].
In the SME, the oscillation probability P νe→νe is given by [35] P νe→νe = |S The transition amplitude S (0) eē is expressed as where U is the PMNS [11] neutrino mixing matrix, E a is the neutrino energy [30], L is the baseline [30], and the sum is over all mass eigenstates a = 1, 2, 3. P (0) is the usual oscillation probability for massive neutrinos in the Lorentz-invariant case. P (1) and P (2) include the interference from common mass-driven mixing and LV mixing. P (1) is calculated as [35] P (1) = 2L · Im [S where (M (1) eē )cd are the experimental factors [30] and δhcd is the LV Hamiltonian. The subscriptcd represents the flavor pairsēē,μμ,ττ ,ēμ,ēτ ,μτ . The experimental factors are defined in terms of the conventional eigenvalues and elements of the PMNS matrix: where τ (1) ab (E, L) =