Study of the wave packet treatment of neutrino oscillation at Daya Bay

The disappearance of reactor $\bar{\nu}_e$ observed by the Daya Bay experiment is examined in the framework of a model in which the neutrino is described by a wave packet with a relative intrinsic momentum dispersion $\sigma_\text{rel}$. Three pairs of nuclear reactors and eight antineutrino detectors, each with good energy resolution, distributed among three experimental halls, supply a high-statistics sample of $\bar{\nu}_e$ acquired at nine different baselines. This provides a unique platform to test the effects which arise from the wave packet treatment of neutrino oscillation. The modified survival probability formula was used to fit Daya Bay data, providing the first experimental limits: $2.38 \cdot 10^{-17}<\sigma_{\rm rel}<0.23$. Treating the dimensions of the reactor cores and detectors as constraints, the limits are improved: $10^{-14} \lesssim \sigma_{\rm rel}<0.23$, and an upper limit of $\sigma_{\rm rel}<0.20$ is obtained. All limits correspond to a 95\% C.L. Furthermore, the effect due to the wave packet nature of neutrino oscillation is found to be insignificant for reactor antineutrinos detected by the Daya Bay experiment thus ensuring an unbiased measurement of the oscillation parameters $\sin^22\theta_{13}$ and $\Delta m^2_{32}$ within the plane wave model.


Neutrino oscillation in the plane wave approximation
The neutrino, a light electrically neutral fermion participating in weak interactions, was suggested by Pauli to save the conservation of energy and momentum in nuclear β-decays. Since then, three flavors of neutrinos ν α = (ν e , ν µ , ν τ ) were discovered, each produced or detected in association with a corresponding lepton α = (e, µ, τ ). The neutrinos, which are completely parity-violating in their weak interactions, suggested that the gauge group of the electro-weak sector of the remarkably successful Standard Model (SM) should be built using fermions with left-handed chirality. Given the unique properties of neutrinos, studies of them may reveal a path to physics beyond the SM. In the past, experiments observing solar and atmospheric neutrinos brought increased attention to neutrino physics due to long-standing discrepancies between detection rates and no-oscillation models. Despite an impressive number of proposed solutions to these problems, all were successfully resolved by the hypothesis of neutrino oscillation, first proposed by Pontecorvo [1,2] in the late 1950's. Neutrino oscillation is a phenomenon firmly established in experiment, which has been observed with solar [3,4,5], atmospheric [6,7], particle accelerator [8,7] and reactor [9,10,11,12] neutrinos.
Neutrino oscillation is a quantum phenomenon of quasi-periodic change of neutrino flavor ν α → ν β with time. This phenomenon originates in the nonequivalence of neutrino flavor ν α and mass ν k = (ν 1 , ν 2 , ν 3 ) eigenstates, differences in their masses, and an assumption that the produced and detected neutrino states are coherent superpositions of neutrino mass eigenstates: where V αk is an element of the unitary PMNS-matrix, named after Pontecorvo, Maki, Nakagawa, Sakata, and p is the momentum of the neutrino. The time evolution of the state in Eq. (1) is expressed as where E k = p 2 + m 2 k . This leads to the oscillatory behavior of the probability to detect a neutrino originally of flavor α as having flavor β: where L osc kj = 4πp/∆m 2 kj is the oscillation length due to the non-zero differences ∆m 2 kj = m 2 k − m 2 j , and time t is approximated by the traveled distance L. The underlying theory, assuming a plane wave approximation, was developed in the middle of the 1970s [13,14,15]. Although successful in explaining a wide range of neutrino experiments, it is well known that this approximation is not self-consistent, and leads to a number of paradoxes [16,17]. The applicability of the plane wave approximation is discussed in detail in Refs. [18,16,19,20]. After the first theory was developed, Refs. [21,22,23,24] pointed out the necessity of a wave packet treatment of neutrino oscillation.

Wave packet treatment of neutrino oscillation
The wave packet is a coherent superposition of different waves whose momenta are distributed around the most probable value, with a certain "width" or dispersion. Therefore, a wave packet is localized in space-time as well as in energy-momentum space. The wave packet formalism facilitates the resolution of the paradoxes of the plane wave theory, and predicts the existence of a coherence length. The latter arises due to the different group velocities of a pair ν k and ν j , which causes a separation in space over time. The smallness of the differences of neutrino masses relative to their typical energies suggests that the coherence length of neutrino oscillation is the largest among all known phenomena.
After the pioneering studies [21,22,23], the wave packet models of neutrino oscillation were developed in roughly two varieties. The first one relies on a relativistic quantum mechanical (QM) formalism that does not predict the dispersion of the neutrino wave packet in momentum space, such as in Refs. [18,19,25]. The second one is based on calculations within quantum field theory (QFT), describing all external particles involved in neutrino production and detection as wave packets while treating neutrinos as virtual particles. The neutrino wave-function is then calculated rather than postulated. The effective momentum dispersion of the neutrino wave function depends on the kinematics of neutrino production and detection and on the momentum dispersions of the external particles, as in Refs. [26,27,28,29,30,31,32]. Both approaches predict a number of observable effects, like a quantitative condition on the coherence of mass eigenstates in the production-detection processes, as well as a loss of coherence.
In wave packet models, the intrinsic momentum dispersion σ p of the neutrino wave packet is an effective quantity comprising the microscopic momenta dispersions of all particles involved in the production and detection of the neutrino. A non-zero value of σ p leads with time to the decoherence in the quantum superposition of massive neutrinos which results in a vanishing oscillation pattern of ν α → ν β transitions. In addition, the oscillation pattern is smeared further in the reconstructed energy spectrum due to a non-zero experimental resolution δ E of the neutrino energy.
Despite considerable progress in building wave packet models, none of these approaches provides a solid quantitative theoretical estimate of σ p or of the spatial width σ x = 1/2σ p . Theoretical estimates vary by orders of magnitude, associating the dispersion of the neutrino wave packet with various scales; for example, uranium nucleus size (σ x 10 −11 cm, σ p 1 MeV), atomic or interatomic size (σ x (10 −8 − 10 −7 ) cm, σ p (10 3 − 10 2 ) eV), pressure broadening (σ x 10 −4 cm, σ p 0.1 eV), etc. While the current literature does not include calculations of the neutrino wave function from first principles for any type of neutrino experiment 2 , it also lacks experimental investigations of decoherence effects in neutrino oscillation inferred from the finite size of the neutrino wave function 3 .
One of the motivations of this paper is to provide a first search for a possible loss of coherence in the quantum state of neutrinos following from the wave packet treatment of neutrino oscillations, using data from the Daya Bay Reactor Neutrino Experiment. The second motivation is to demonstrate that the oscillation parameters estimated with the plane wave approximation are unbiased. The oscillation probability formula modified by the wave packet contribution, which is discussed further, has two distinctive features: it depends on ∆m 2 kj /p 2 σ rel via the so-called localization term and on L∆m 2 kj σ rel /p via the term responsible for the loss of coherence with distance. The large statistics, good energy resolution, and multiple baselines of the Daya Bay experiment make its data valuable in the study of these quantum decoherence effects in neutrino oscillation.

Neutrino oscillation in a wave packet model
Measured energy spectra ofν e interactions are compared to a prediction using a QM wave packet model of neutrino oscillation which is briefly outlined in what follows. We simplify the consideration by examining a one-dimensional wave packet of the neutrino 4 . The plane wave state in Eq. (1) is replaced by a wave packet describing a neutrino produced as flavor α: is the wave function of the neutrino in momentum space and is assumed to be Gaussian: where the subscript P in f P (p), p P and σ pP indicates the quantities at production. In configuration space the state in Eq. (4) describes a wave packet with mean coordinate x P at time t P . The state in Eq. (4) is normalized as ν α (p P ; t P , x P )| ν α (p P ; t P , x P ) = 1. Similarly, a wave packet state at detection | ν β (p D ; t D , x D ) is defined as the state given by Eq. (4).
A projection of | ν α (p P ; t P , x P ) onto ν β (p D ; t D , x D )| produces the flavorchanging amplitude which depends on L ≡ x D −x P , time difference t D −t P and on the effective mean neutrino momentum p and momentum dispersion σ p comprising the details of production and detection 5 The probability |A αβ (p; t D − t P , L, σ p )| 2 should be integrated over production time t P (or, equivalently, over t D − t P ) and most probable momentum p P to get an experimentally observable oscillation probability: where the phase ϕ kj is the sum of the plane wave phase ϕ kj = 2πL/L osc kj and correction ϕ d kj due to the dispersion of the wave packet: Oscillation probability formulas similar to Eq. (8) but neglecting wave packet dispersion were obtained in several studies (see, for example, Refs. [29,18,31,39]). Eq. (8) has appeared as a particular case of a more general consideration within QFT with relativistic wave packets [32]. Relativistic invariance suggests that σ rel should be Lorentz invariant. In the QM approach adopted in Eq. (4)-Eq. (8) the only possibility to preserve Lorentz invariance is for σ rel to be a constant 6 . The probability in Eq. (8) contains three quantities with dimensions of length: where σ rel = σ p /p, L osc kj is the usual oscillation length of a pair of neutrino states |ν k and |ν j , L coh kj is interpreted as the neutrino coherence length, i.e. the distance at which the interference of neutrino mass eigenstates vanishes, and finally L d kj is the dispersion length, i.e. a distance at which the wave packet is doubled in its spatial dimension due to the dispersion of waves moving with different velocities. The term suppresses the coherence of massive neutrino states |ν k and |ν j if ∆m 2 kj σ m 2 , where σ m 2 = 2 √ 2pσ p could be interpreted as an uncertainty in the neutrino 6 Since the QFT approach considers both neutrino production and detection one finds that σ rel , being a relativistic invariant, is actually a function of kinematic variables involved in the production and detection processes as well as of momentum dispersions of wave packets describing all involved particles [40]. Therefore, in comparing the QM and QFT approaches, we may treat the QM σ rel as that of the QFT approach averaged over the kinematic variables of all external wave packets involved in neutrino production and detection. mass squared [22]. D 2 kj can be seen from another perspective as the localization term suppressing the oscillation if √ 2πσ x L osc kj , where σ x = (2σ p ) −1 is the width of neutrino wave packet in the configuration space.
It is worth mentioning that terms in Eq. (8) which correspond to the interference of ν k and ν j states also get suppressed by the denominator and vanish for both limits σ p → 0 and σ p → ∞, reducing the oscillation probability in Eq. (8) to the non-coherent sum which does not depend on energy and distance.
For theν e at Daya Bay, 1 − P ee is expressed as

Sensitivity of Daya Bay experiment to neutrino wave packet
The Daya Bay experiment is composed of two near underground experimental halls (EH1 and EH2) and one far underground hall (EH3). Each of the experimental halls hosts identically designed antineutrino detectors (ADs). EH1 and EH2 contain two ADs each, while EH3 contains four ADs. Electron antineutrinos are produced in three pairs of nuclear reactors via β decays of neutron-rich daughters of the fission isotopes 235 U, 238 U, 239 Pu and 241 Pu, and detected via the inverse β decay (IBD). The coincidence of the prompt (e + ionization and annihilation) and delayed (n capture on Gd) signals efficiently suppresses the backgrounds, which amounted to less than 2% (5%) of the IBD candidates in the near (far) halls [41]. The Gd-doped liquid scintillator target is a cylinder of three meters in both height and diameter. The detectors have a light yield of about 165 photoelectrons/MeV and a reconstructed energy resolution δ E /E ≈ 8% at 1 MeV of deposited energy in the scintillator. More details on the experimental setup are contained in Refs. [41,42,43,44].
The studies in this paper are based on data acquired in the 6-AD period when there were two ADs in EH1, one AD in EH2 and 3 ADs in EH3, with the addition of the 8-AD period from October 2012 to November 2013, a total of 621 days. The number of IBD candidates used in this analysis, and the mean baselines of the three experimental halls to each pair of reactor cores, are summarized in Table 1  of the reactor-to-target expectation with the detector-response function. The reactor-to-target expectation takes into account the antineutrino fluxes from each reactor core including non-equilibrium and spent nuclear fuel corrections, first order in 1/m p (m p =proton mass) IBD cross-section accounting for the positron emission angle [45], and the oscillation survival probability P ee given by Eq. (3) for the plane wave model and by Eq. (8) for the wave packet model. The detector response-function accounts for energy loss in the inner acrylic vessel, liquid scintillator and electronics non-linearity and energy resolution δ E . For relatively large values of σ p δ E , the effects of these two parameters on the observed energy spectra might appear similar, however they are distinct. First, they have different physical origins: while σ p is governed by the most localized particle in the production and detection of the neutrino, δ E is determined by the energy depositions of the final state particles in the detector. Second, these effects can also be distinguished from their order of occurrence since the microscopic processes used in the energy estimation occur later in time with respect to the neutrino interaction in the detector. Third, their effects are not identical. In particular, as described in Sec. 2.1, the limit σ p → 0 leads to the decoherence of neutrino oscillation in contrast to the impact of energy resolution which does not lead to any smearing in the reconstructed energy spectrum in the limit δ E → 0.
In order to illustrate analytically an interplay of σ p and δ E , let us consider the exponential in the oscillation probability in Eq. (8) convolved with a Gaussian energy resolution, as a function of the reconstructed energy E vis , assuming δ E p, infinite dispersion length L d , neglecting the D 2 term, and suppressing mass eigenstate indices for the sake of compactness 7 : where L osc and L coh are given by Eq. (10) and the effective coherence length comprises both the intrinsic σ p and detector resolution δ E : where L osc rec and L coh rec are given by L osc and L coh replacing p with E vis , and L coh det is given by L coh rec , replacing σ p with δ E . The interplay of σ p and δ E is illustrated by the effective coherence length L coh eff , which is dominantly determined by the smallest among L coh rec and L coh det , or by the largest among σ p and δ E . For illustrative purposes Fig. 1 shows the ratio of the observed to expected numbers of IBD events assuming no oscillation using the data collected at the near and far experimental halls as a function of reconstructed visible energy E vis . Figure 1 also shows the expected ratio for neutrino oscillation with the plane wave and wave packet models with σ rel of 0.33 and 8 · 10 −17 as examples.
Both model expectations are shown with the oscillation parameters fixed to their best-fit values within the plane wave model 8 . For this set of parameters, the wave packet models with σ rel = 0.33 and with σ rel = 8·10 −17 are inconsistent 7 The actual implementation of the detector effects in this analysis was performed numerically without approximations 8 The following values of the oscillation parameters were used in Fig. 1: ∆m 2 21 = 7.53 · 10 −5 eV 2 , ∆m 2 32 = 2.45 · 10 −3 eV 2 , sin 2 2θ 12 = 0.846, sin 2 2θ 13 = 0.0852.
with the data by about five standard deviations, thus motivating the chosen values of σ rel . The two panels illustrate how the visible energy spectra are modified in the near and far halls depending on the intrinsic dispersion of the neutrino wave packet. Remarkably, most changes in the energy spectra due to σ rel are in opposite directions for near and far halls, which can be explained qualitatively as follows. As mentioned above, the extremes σ p → 0 and σ p → ∞ would yield fully decoherent neutrinos with the oscillation probability given by Eq. (12). Antineutrinos detected at the near halls experience a relatively small oscillation in the plane wave approach. The values of σ rel selected for Fig. 1 make theν e partially decoherent and P ee tend towards Eq. (12), predicting a smaller number of survivingν e as compared to the plane wave formula. The distance at which the far detectors of the Daya Bay experiment are placed is tuned to observe the maximal oscillation effect due to ∆m 2 32 . Partial decoherence of theν e tends to reduce the oscillation, thus predicting a larger number of survivedν e with respect to the plane wave formula. This feature of Daya Bay provides additional sensitivity to the decoherence effects and makes such a study less sensitive to the predicted reactorν e spectrum.
These results demonstrate that one could obtain reasonable fits of the data within the wave packet model with certain values of σ rel and yield best-fit values of the oscillation parameters which differ from the corresponding best-fit values with the plane wave model, assuming normal mass hierarchy 9 : ∆m 2 32 = 2.45 · 10 −3 eV 2 , sin 2 2θ 13 = 0.0852, χ 2 /ndf = 245.9/(256 − 3).
However, Eqs. 16, 17 do not correspond to the global minimum of the χ 2 discussed below because σ rel was fixed to two arbitrary values for illustrative purposes. In order to find the global minimum we performed a detailed statistical analysis of the allowed region of σ rel .

Statistical framework
As the goodness-of-fit measure we use , where d is a data vector containing detected numbers of IBD candidates in energy bins and in different detectors, while t(η) is the corresponding theoretical model vector which depends on constrained and unconstrained parameters η.
All constraints of the model as well as expected fluctuations in the number of IBD events are encompassed in the covariance matrix V . The model vector t(η) comprises expected numbers of IBD and background events. All constrained parameters (or systematic uncertainties) relevant for the Daya Bay oscillation analyses were taken into account in this analysis. These are mainly associated with the reactor antineutrino flux, background predictions and the detector response modeling. The uncertainty of the detector response is dominant. Details can be found in Refs. [41,44].
The analysis was done with four unconstrained parameters σ rel , ∆m 2 32 , sin 2 2θ 13 and reactor flux normalization N . The confidence regions are produced by means of two statistical methods: the conventional fixed-level ∆χ 2 analysis and the Feldman-Cousins method [46]. The marginalized ∆χ 2 statistic is where η = (σ rel , ∆m 2 32 , sin 2 2θ 13 , N ) and η is its subspace with parameters of interest (η = σ rel for one dimensional interval, and η = (σ rel , ∆m 2 32 ) or η = (σ rel , sin 2 2θ 13 ) for two dimensional regions), and both are used to determine the p-value of the observed dataset and the model.
The closed interval corresponding to the 100·(1−α)% confidence level (C.L.) is constructed for both the fixed-level ∆χ 2 analysis and the Feldman-Cousins method as the region of η which satisfies: where ∆χ 2 1−α is the (1−α)-th quantile of the statistic in Eq. (19). The tabulated values of the quantile χ 2 n;1−α of the χ 2 n distribution with n degrees of freedom (n = 1, 2 for one and two dimensional confidence regions) were used for the fixed-level ∆χ 2 analysis. Toy Monte Carlo sampling was used to determine ∆χ 2 1−α of the statistic in Eq. (19) with the Feldman-Cousins method. An open confidence interval can be constructed if neutrinos are assumed to be produced and detected coherently, which is equivalent to assuming σ rel 10 −16 . In this case, instead of using Eq. (19), an upper bound on σ rel can be computed using the modified statistic [47] withσ rel representing the best-fit value. In the fixed-level ∆χ 2 analysis the 100 · (1 − α)% C.L. upper limit is given by: For example, in order to set a 95% C.L. upper limit, the quantile χ 2 1;0.9 = 2.71 was used. The Feldman-Cousins method automatically produces the proper interval using the interval construction in Eq. (20). Figure 2 displays the allowed regions in (∆m 2 32 , σ rel ) and (sin 2 2θ 13 , σ rel ) obtained with both the fixed-level ∆χ 2 and the Feldman-Cousins methods, which are found to be consistent. For the values of σ rel 10 −16 the decoherence effects lead to strong correlations between ∆m 2 32 , sin 2 2θ 13 and σ rel , yielding smaller values of ∆m 2 32 and larger values of sin 2 2θ 13 . These correlations are expected taking into account the explicit form of 1 − P ee (L) in Eq. (13). For σ rel O(0.1), these correlations are found to be significantly weaker.
The upper bound of Eq. (24) corresponds to L coh 32 > 1.94 L osc 32 /2 and L d 32 > 2.96 L osc 32 /2. The lower bound can also be interpreted in terms of length σ x which corresponds to the spatial width of the neutrino wave packet. Taking the average momentum p = 4 MeV of detected reactorν e , the lower bound of Eq. (24) rules out σ x 1 km. The Daya Bay data is not sensitive enough to constrain the D 2 kj term significantly better. Thus, the lower limit is much weaker than an obvious constraint of σ x 2 m which follows from the consideration that the σ x (which equals 1/2σ p ) ofν e wave packets detected by the Daya Bay Experiment does not exceed the dimensions of the reactor cores and detectors. Taking this constraint into account, σ p 5·10 −8 eV, which for the average momentum p = 4 MeV, translates into σ rel 10 −14 . Such a σ rel corresponds to the regime where D 2 kj 1 and the localization term can be safely neglected, thus allowing us to put an upper limit of: σ rel < 0.20, at a 95% C.L.

Summary
We performed a search for the footprint of the neutrino wave packet which should show itself through specific modifications of the neutrino oscillation probability. The reported analysis of the Daya Bay data provides, for the first time, an allowed interval of the intrinsic relative dispersion of neutrino momentum 2.38 · 10 −17 < σ rel < 0.23. Taking into account the actual dimensions of the reactor cores and detectors, we find that the lower limit σ rel > 10 −14 corresponds to the regime when the localization term is vanishing, thus allowing us to put an upper limit: σ rel < 0.20 at a 95% C.L. The obtained limits can be read as 10 −11 cm σ x 2 m.
The current limits are dominated by statistics. With three years of additional data the upper limit on σ rel is expected to be improved by about 30%. The allowed decoherence effect due to the wave packet nature of neutrino oscillation is found to be insignificant for reactor antineutrinos detected by the Daya Bay experiment thus ensuring an unbiased measurement of the oscillation parameters sin 2 2θ 13 and ∆m 2 32 within the plane wave model. WP (σ rel = 0.33) Figure 1: Ratios of the observed to expected numbers of IBD events in the absence of oscillation as a function of reconstructed visible energy E vis . The data are grouped by near (EH1+EH2) and far (EH3) halls, displayed in the upper and in the bottom panels respectively, with the error bars representing the statistical uncertainties. Superimposed solid lines are ratios assuming neutrino oscillations within the plane wave model (PW) with the best-fit values of sin 2 2θ 13 and ∆m 2 32 obtained with the plane wave model. The ratios using the wavepacket model (WP) assume σ rel = 0.33 (dashed line) and σ rel = 8 · 10 −17 (dot-dashed line), as two examples. The green lines correspond to the wave packet model ratios assuming the bestfit values of sin 2 2θ 13 and ∆m 2 32 obtained with the plane wave model and thus, inconsistent with the data by about five standard deviations. The red lines correspond to the wave packet model ratios assuming the best-fit values of sin 2 2θ 13 and ∆m 2 32 obtained within the wave packet model, yielding a much better agreement with the data. All ratios enter the region below 2me, which corresponds to the IBD threshold, because of detector response effects like energy reconstruction and absorption in the inner acrylic vessel (see details in Refs. [41,44]).  : Allowed regions of (∆m 2 32 , σ rel ) (top) and of (sin 2 2θ 13 , σ rel ) (middle) parameters obtained with fixed-level ∆χ 2 (contours corresponding to 1σ, 2σ, 3σ C.L., dashed lines) and within the Feldman-Cousins (contours corresponding to 1σ, 2σ C.L., solid lines) methods. Bottom panel shows the marginalized ∆χ 2 (σ rel ) statistic given by Eq. (19) vs σ rel . Note the break in the abscissa and the change from a logarithmic to linear scale.