Neutrino Oscillation Measurements Computed in Quantum Field Theory

We perform a calculation in quantum field theory of neutrino oscillation probabilities, where we include simultaneously the source, detector, and neutrino fields in the Hamiltonian. Within the appropriate limits associated with current neutrino oscillation experiments, we recover the standard oscillation formula. On the other hand, we find that the dominant contributions to the amplitude are associated with different neutrino mass eigenstates being emitted at different times, such that they arrive at the detector at the same time. This is contrary to the neutrino wave packet picture, where they are emitted simultaneously and separate as they travel to the detector. This has direct consequences regarding the mechanisms that lead to a damping of neutrino oscillations for very long baselines. Our analysis also provides a pedagogical example of a measurement process in quantum mechanics.


I. INTRODUCTION
The experimental observation of flavor transitions in neutrino oscillation experiments implies that neutrinos have mass. The neutrino field ν which couples to the charged lepton is, in general, a linear combination of the fields ν a which diagonalize the neutrino mass matrix, Global combinations of data from oscillation experiments find good agreement using the following formula for the probability of neutrino oscillation [1,2] from flavor to : where U a is the mixing matrix in Eq. (1), m a is the mass of the neutrino mass eigenstate ν a , L is the baseline of the experiment (i.e., the distance between source and detector) and E ν is the neutrino energy.
Here, ξ is a constant, which is typically derived to be ξ = 1/2, and the neutrino oscillation parameters obtained from oscillation experiments are quoted using this value in Eq. (2). An overview of the standard derivation of the neutrino oscillation phase that leads to ξ = 1/2, e.g., the one discussed in the PDG [3], can be found in Appendix A. Here, the oscillation probability P(T, L) is not only a function of the baseline distance L, but also of the transit time T that neutrinos take between their production and detection. It is then imposed, using an ultra-relativistic limit, that T L for all neutrino species, and the result is ξ = 1/2. This discussion can be extended to allow for the damping of oscillations for very long baselines L. Many discussions exist in the literature regarding this phenomenon (for example, see Refs. ). Common to these treatments is that all neutrino mass eigenstates are produced by the source particle in the same region of time. Under this condition, there is an unavoidable source of damping of oscillations for large enough values of L, since, due to their different group velocities, the different mass components of the neutrino wave packet cease to overlap as they undergo their journey from source to detector. The resulting damping of neutrino oscillations is estimated to be negligible for terrestrial experiments.
There is an inconsistency, however, when computing neutrino oscillations in this manner. The neutrino transit time T between production and detection is an unmeasured variable in contemporary oscillation experiments. One can only calculate the probability for the outcome of performed experiments. If T is unmeasured, one cannot calculate the outcome for such an experiment as a function of T . Instead, all amplitudes with different values of T must be summed. We can represent this with the following cartoon: Here, the total amplitude A is the sum over all unmeasured variables, namely the time t that the source particle (depicted by a purple line) decays to a neutrino (blue lines) and other particles (green lines), and the mass of the neutrino m a . We have drawn two mass eigenstates, i.e., a = 1, 2. The vertical direction is time, and the horizontal direction is space. The yellow dot represents the spacetime region of a detector. As a result of our calculation, we find that the following diagrams are the dominant contribution to the amplitude after integrating over time, i.e., they are the extrema of the path integral: The oscillation probability P is the square of this amplitude, P = |A| 2 . Illustrating the result when one calculates this sum over amplitudes is the primary purpose of this paper. Given the relevant approximations for typical neutrino oscillation experiments, we find that the amplitude is dominated by contributions where the source particle decays at different times such that all the different neutrino species interfere with the detector at the same time.
In this paper, neutrino oscillation experiments will be treated entirely within quantum field theory, including simultaneously the source, neutrinos, and detector. To simplify the algebra, the source and detector particles are treated using heavy particle effective theory, where, to lowest order in 1/M , the particles can be treated as Wilson lines with a constant velocity v µ . The results we obtain in this paper hold even without this assumption, but the calculations are more involved (they can be found in Appendix B). The main difference from previous analyses is a full quantum mechanical treatment of the problem. We start with the source in an excited state and the detector in the ground state at t = 0, and compute the probability to find the detector in an excited state at a later time. The source and detector make transitions between their excited and ground states by emitting and absorbing neutrinos. We make no assumptions about the neutrino, or about its journey from source to detector. In quantum mechanics, one cannot make assumptions about the intermediate state, since it is not measured. In the Feynman path integral approach, the neutrino takes all possible journeys between source and detector which have to be summed over, and the different mass eigenstates do not have to travel together in the path integral.
We work in the limits E ν m a , E ν 1/L, where the detection amplitude can be estimated analytically. The amplitude is dominated by a semiclassical trajectory, which is determined by the structure A schematic diagram of the dominant contributions to the amplitude for neutrino propagation at oscillation experiments. We show two mass eigenstates ν 1 (light blue) and ν 2 (dark blue), where ν 1 is the lighter neutrino. The light and dark green lines correspond to other decay products produced when the source decays to the neutrino. The ν 1 and ν 2 contributions to the amplitude have to be added before squaring the amplitude to get the probability. The source (purple lines) is stationary in the figure, as occurs when neutrinos are produced in reactor experiments. This configuration is the sum of dominant amplitudes when ν 1 and ν 2 are emitted from the source at different times, so that they arrive at the detector at the same time.
of the amplitude integral. It is important to keep in mind that this trajectory depends on the entire experiment, i.e., on the properties of the source, detector, and neutrino, and should not be taken too literally. If the source amplitude evolves in time as a pure exponential e −i∆ S t−Γ S t/2 , we find that the dominant semiclassical trajectory is the one shown in Fig. 1. The two neutrino mass eigenstates ν 1,2 are emitted at different times from the source, in such a way that they arrive at the detector at the same time. Figure 1 cannot be taken literally as a spacetime diagram of the neutrino oscillation experiment. It represents the semiclassical trajectory which dominates the amplitude for the calculation that we are doing in this paper. A different calculation with a different time-dependence for the source, such as a source which is turned on and off, requires recomputing the large L amplitude, and the result is not determined by using Fig. 1.
The outline of this paper is as follows. In Sec. II, we discuss the HQET-based formalism we use for neutrino sources. Section III computes the oscillation probability for a stationary source and detector treated as two-state systems. Section IV generalizes the calculation to a moving source and stationary detector. We show how the picture in Fig. 1 follows from the calculations in Secs. III, IV, and discuss the implications of our results in Sec. V. Appendix A summarizes the factor of two puzzle for ξ. Appendix B shows how our results can be generalized, by including localization effects and recoil corrections for the source, and how the results can be applied to sources other than two-state systems. In particular, the results of the appendices show that calculations are applicable for both reactor neutrino experiments, and accelerator experiments where the neutrinos are produced by pion beams.

II. HEAVY SOURCES AND DETECTORS OF NEUTRINOS
In oscillation experiments, neutrinos are produced by the decay of a quasi-stable source particle, which is sufficiently localized within a region of space. 1 If there exists a reference frame in which the source is non-relativistic, and remains so through time scales of order its lifetime, then such a state can systematically be treated within heavy particle effective field theory. To illustrate, the amplitude for a free, non-relativistic, scalar particle with mass M to propagate from position eigenstate x to x is: Here, τ is the time in the particle's rest frame, and H 0 M + p 2 /2M . The limit M → ∞, i.e., ignoring additional terms in Eq. (3), neglects the spread of the source's spatial wave function in time due to free evolution. Therefore, the state can be taken to be arbitrarily localized in the M → ∞ limit. A useful mnemonic is: so that heavy particles in quantum mechanics can simultaneously have definite position and velocity, but not definite position and momentum [28]. For the purposes of our present analysis, we will be treating the source as a particle with ∆x ∼ 0 and ∆v ∼ 0. Corrections can be included in the 1/M expansion, just as in heavy particle effective theory, but such corrections are negligible for terrestrial neutrino experiments. 2 We now explain the Hamiltonian we will use to compute the neutrino oscillation probability. The system consists of a source S, a detector D and neutrino fields. The source can transition to a neutrino ν and a collection of other particles X. In the heavy particle limit, the source travels along a straight worldline [29][30][31] from which neutrinos and other particles are produced (see Appendix B for more details). The source S behaves like a B meson in HQET [31]. The S → νX decay can be treated like B → D decays if X is also a heavy state, as in nuclear β decay, or like B → π decays if X contains light states, as in π → ν decay. In either case, we can still use a 1/M S expansion for the source, which simplifies the integrals that need to be evaluated.

III. OSCILLATIONS WITH A STATIONARY SOURCE
Here, we consider a scalar theory with a (heavy) stationary source, located at the origin, which can have two quantum states. The source is originally in the excited state S E , which decays to the ground 1 If this were not the case, e.g., if the state of the source particle were well approximated instead by a momentum eigenstate, its spatial wave function would be approximately a constant over all space, and it would not be possible to observe oscillations as a function of the experiment's baseline distance [19]. state S G with no recoil, and creates a quantum of the scalar field φ a with mass m a . 3 ∆ S is the mass splitting of the two source states. Since the source can decay from S E → S G +φ, interactions of the source with the φ field produce a width for the source in the excited state. The well-known Weisskopf-Wigner result [32] shows that this interaction can be included by the replacement ∆ S → ∆ S − iΓ S /2, where Γ S is the decay rate of the source [18,26,33]. We treat the detector in a similar fashion, where the heavy detector particle (called D G for ground state), located at position L, can absorb a field φ a , and transition to another heavy state (called D E for excited state) with no recoil, where ∆ D is the detector mass splitting. This model is similar to the one used in Ref. [34], and is essentially the famous Fermi two-atom problem [35]. See endnote 4 of Ref. [34] for the history of this problem. In order to model a detector that can measure more than one energy, many distinguishable two-level systems with different energy splittings ∆ D j can be placed at the same location. The detector width Γ D can be included in the same way as the source width Γ S . We drop it, since it does not lead to any new features in the calculations.
In Appendix B 3, we discuss the general calculation, in the full quantum field theory, which leads to the same final result as the one presented in this section, but with considerable increase of complexity.
There are several scalar fields φ a with different masses m a . We consider the simple case where the source and detector couple to one linear combination of mass eigenstates φ = a U a φ a , with a |U a | 2 = 1. It is trivial to generalize to the case where source and detector couple to different linear combinations, by a replacement of the mixing matrix U a with ones for the source and detector. The interaction Hamiltonian is where, and E a p ≡ p 2 + m 2 a . The source and detector coupling are g. We consider sources and detectors which are active for a finite time interval, by including the source and/or detector terms in Eq. (5) only during the active time-window. φ (+) annihilates φ, and φ (−) creates φ. The use of only φ (±) in the two terms of V is referred to as the rotating wave approximation. Deviations from this approximation are suppressed due to energy conservation.
The initial quantum state of the system at t = t i = 0 is |i = |S E ⊗ |D G ⊗ |0 φ -the source is in the excited state, the detector is in the ground state, and there are no φ quanta excited. We can set up the initial conditions |i , since we have an infinite amount of time before starting the experiment. The initial density matrix is thus To study whether neutrinos have oscillated, we evolve ρ i in time using Eq. (5), and trace with the final state density matrix at time t f This gives the probability that the detector is in an excited state. It places no restriction on the states of the source and φ field, which are summed over. In our example, the source and detector are on for long times, and the neutrino time of flight is also large compared to the instrinsic time scales 1/∆ S,D . Approximate energy conservation implies that the only source state which contributes is |S G , and the amplitude of |S E is exponentially small. With our Hamiltonian, the only φ state which contributes is the ground state |0 φ . 4 Thus we will study the transitions If there are several detectors, then all detectors are in the ground state in |i , and we have a set of final states |f j where detector j is excited, and all other detectors are in the ground state. Final state configurations where multiple detectors are excited are higher order in the weak interactions, and can be neglected.
The lowest-order contribution in perturbation theory in g that gives rise to the transition between these initial and final states is where the source is on for all time, and the detector is on only during the time window [τ, τ + T ]. τ is the time of detection, and T is the size of the time bin. This allows us to model an experiment where detector events are counted as a function of time. A set of such detectors, which are on at different times [τ 1 , τ 1 + T ], [τ 2 , τ 2 + T ], etc., allows one to compute the detection rate as a function of arrival time at the detector. To maintain a simpler notation in the intermediate steps, we calculate the transition amplitude in Eq. (11) with ∆ S as the energy splitting of the source, and use a single φ mass eigenstate. The full answer is then given by the replacement ∆ S → ∆ S − iΓ S /2 and summing over masses at the end. Here, Eq. (11) becomes where Performing the t 2 and angular integrals, and changing integration variables to E p : The contours used to estimate the value of the energy integral in Eq. (14). The vertical parts of the contour are negligible in the limit ∆ S 1/L, ∆ S m, and ∆ S 1/|t 1 − L|.
The integrand has no poles along the real E p axis for E p > m, but the individual terms do. We therefore first deform the contour to go slightly above E p = ∆ S , and then split the integral into four terms given by multiplying out the two factors in square brakets in Eq. (14). We utilize residue methods using the contours in Fig. 2 to compute the energy integral. Depending on the sign of the exponential for large E p , we consider the closed contour in either the upper or lower half plane, such that the exponential is damped along the circular contour at infinity. The original integral is then given by the residue of the pole at ∆ S if the contour is closed in the lower half plane, plus the integral along the (dotted) vertical section, or by the (solid) vertical section if the contour is closed in the upper half plane. The vertical parts of the contours, parallel to the imaginary axis, are suppressed by powers of 1/L in the limit that ∆ S m, ∆ S 1/L, and ∆ S 1/|t 1 − L|, which is the case in neutrino oscillation experiments. They give a contribution analogous to near-field effects in classical radiation theory.
If L > 0 and t 1 − L < 0, i.e., the detector is outside the source's light cone, the integral vanishes up to the power suppressed terms. There cannot be a signal transmitted by the intermediary φ particle faster than light. In our calculation, there is a small power suppressed contribution outside the light-cone (from the vertical segment) because of the form Eq. (5) where the φ field has been split into its creation and annihilation parts. If we had instead used the full φ field in both terms, causality would be exact, and A j = 0 outside the light-cone (see Ref. [34,36]).
On the other hand, if L > 0 and t 1 − L > 0, and the detector is well inside the source's light cone, the integral is up to power suppressed terms. Making the replacement ∆ S → ∆ S − iΓ S /2, summing over multiple mass eigenstates φ a : The total probability for the transition includes a sum over the distinguishable detector systems labelled by j: where ρ D (∆ D j ) is the number density of detector systems with energy splitting ∆ D j . The sinc function 5 in Eq. (17) should be familiar from the usual quantum mechanical derivation of Fermi's Golden Rule. If Γ S = 0, then in the limit ∆ S T → ∞ as in Fermi's Golden rule. The only detectors which respond are those with the same energy splitting as the source, and the detection probability is proportional to the time T that the detector is on.
For a non-zero width, the detector response is given by a Lorentzian with a linewidth Γ S , and the total probability is not of order T , but of order 1/Γ S , because that is the lifetime of the source. The e Γ S T /2 factor converts e −Γ S (τ +T /2) in |A j | 2 into e −Γ S τ , where τ is the time the detector is first turned on.
In the limit that ∆ S m a , the total probability becomes if Γ S = 0, and is the number of detectors at energy ∆ S , smeared over a width Γ S . Since Γ S ∆ S , Eq. (21) contains the well-known oscillation term, the classical 1/L 2 area law for flux from a point source in three spatial dimensions, as well as an exponential damping factor, which depends on the masses m a . To further simplify the expression, it is useful to note that the velocity of an ultra-relativistic particle in the lab frame with energy E can be approximated as Since ∆ S is the energy of all the intermediary particles φ a , independent of their masses, the probability Eq. (23) has the simple form The prefactor g 2 ρ D /2πΓ S L 2 is associated with the probability of production, and detection, and the inverse square law. The oscillation factor is The derivation of P osc did not assume anything about the journey of the intermediate neutrino state.
In particular we made no mention of any neutrino wavepackets travelling from source to detector. The result Eq. (27) follows by using time-ordered perturbation theory, and keeping the leading terms in the L → ∞ limit.
We can now compare our oscillation result with the standard formula Eq. (2). The mixing angle is |U a | 2 since we have assumed that source and detector both couple to the same linear combination of φ. The imaginary part of the exponential is the same as Eq. (2) with ξ = 1/2. The phase factor is the same as that obtained using p a L, where p a = ∆ 2 S − m 2 a . The damping term is interesting. It says that the detector sees the source at the earlier time τ − L/v a , which depends on the neutrino species a. Thus the different neutrino species are emitted from the source at different times, in such a way that they arrive at the detector at the same time. This is very different from the usual discussion, in which the source emits a neutrino wavepacket at t = 0, which then propagates and spreads so that the different components arrive at the detector at different times. The consequences of Eq. (27) are discussed further in Sec. V. The geometrical picture in Fig. 1 leads to the phase and damping factors in Eq. (27). We derive this in the next section, after we consider a moving source.
In the above calculation, the exponential decay of the source damps the oscillations for baselines beyond where m a and m b are two different neutrino masses. For example, for a pion source producing ∆ S ∼ 1 GeV neutrinos, the ratio of the coherence length L coh to the oscillation length

IV. OSCILLATIONS WITH A MOVING SOURCE
In neutrino oscillation experiments with pion beams, the source moves relativistically in the lab frame. Using the model in Section III, we calculate the transition amplitude with a source which begins at the origin at t = 0, and moves with a velocity v in the lab frame. The detector remains stationary. The relevant terms in the interaction Hamiltonian in Eq. (5) are modified to: The dependence on v and γ ≡ (1 − v 2 ) − 1 2 is determined by special relativity and is derived from a microscopic quantum field theory in Appendix B. The analogous expression for the transition amplitude in Eq. (12) becomes: One can evaluate these integrals following the same procedure outlined in Section III. In order to illustrate a specific example, we take L = Lẑ and v = vẑ. The calculation is a bit tedious in three spatial dimensions due to the interplay between the momentum and angular integrals. In order to illustrate the essential behavior of the oscillation probability, we do the calculation in one spatial dimension: Using the methods of Section III, the amplitude to excite detector j inside the light-cone when L > vτ (so that the source does not hit the detector) is Here E a ≡ γ ∆ S + v ∆ 2 S − m 2 a is the energy of the neutrino in the lab frame. Eq. (31) differs from the one in three spatial dimensions by a factor of ∆ 2 S − m 2 a in the denominator instead of L 2 . There is no fall-off with L in one dimension. Importantly, because the energy of the neutrino now depends on its mass, our ideal model of a detector is capable, in principle, of distinguishing these different energies, and oscillations need not occur. For this to be possible, we need Γ S and 1/T to both be much smaller than the neutrino energy difference. Real-world detectors are far from ideal, so if there are two neutrinos with energy difference of order (m 2 2 − m 2 1 )vγ/∆ S , and this difference is small enough to excite the same detector particle, then oscillations can occur. The oscillation probability can be estimated as in the previous section, making the replacement ∆ S → ∆ S − iΓ S /2 expanding in m 2 a , using Γ S ∆ S , and picking out the term analogous to Eq. (27) using Eq. (20), where we have assumed that the neutrino energy difference is small compared with Γ S , i.e. the experimental setup cannot resolve the energy difference. Equation (32) reduces to Eq. (2) with ξ = 1/2 in the limit that τ ≈ L and Γ S /∆ S is small, with ∆ S replaced by the relativistic Doppler shifted value We can set τ ≈ L because if τ L, the damping factor makes the amplitude negligible. The exponential decay factor in Eq. (32) has a simple geometrical interpretation, as in Fig. 1. Γ S /γ is the inverse lifetime of the source in the lab frame. Using straight line paths for the source and neutrino with velocity v and v a respectively, the emission time t a (i.e. the time component of x 1 or x 2 in the figure) is The neutrino velocity is Substituting Eq. (35) in Eq. (34) and expanding in m 2 a gives the coefficient of Γ S /(2γ) in Eq. (32). The neutrino propagation phase from Fig. 1 is Expanding this in m a gives which is twice the phase in Eq. (32). However, there is another contribution to the neutrino phase from the phase of the source The first term is independent of neutrino type. The second term cancels half the contribution in Eq. (37), to give the standard oscillation result in Eq. (2) with ξ = 1/2. Eq. (37) is the ξ = 1 calculation outlined in Appendix A, and is the correct neutrino propagation phase. However, the total neutrino phase at the detector also involves the phase of the source, which converts ξ = 1 to ξ = 1/2. The usual derivation in the literature [3] of the ξ = 1/2 value outlined in Appendix A ignores the phase of the source, and assumes neutrinos are emitted at the same time.

V. DISCUSSION AND CONCLUSIONS
The source and detector in neutrino oscillation experiments can be systematically treated using heavy particle effective field theory. In the heavy-particle limit, the source and the detector follow straight paths through spacetime, along which neutrinos can be created or absorbed. Corrections to this treatment of the source can systematically be applied using the 1/M expansion, and in general can lead to a damping of the oscillations. However, the baseline L required to observe such damping from these corrections is many orders of magnitude larger than the diameter of Earth, so this is a negligible effect for contemporary neutrino oscillation experiments.
Oscillation experiments are typically performed with neutrinos that propagate a large spatial distance L such that E ν 1/L. The neutrino propagator is then dominated by its pole. This allows us to illustrate the salient features of neutrino oscillation, treating the neutrino as a scalar and the source and detector system as well-localized, two-level systems, which undergo no recoil upon emitting or absorbing a neutrino. Including the source, neutrino, and detector simultaneously, and performing standard timedependent perturbation theory, we calculate the oscillation behavior both when the source particle is stationary in the lab frame and when it moves. In the limit that E ν m ν , E ν 1/L, and the source does not traverse a large distance (compared to the oscillation length scale) before it decays, we recover the standard oscillation formula used by oscillation experiments. Our treatment of neutrino oscillations can also be applied to mesonic oscillations. Our calculations show that the detector sees the source at different earlier times, due to the different masses of the neutrinos. This differs from statements frequently found in the literature, e.g., Refs. , where it is said that all neutrino mass eigenstates are produced at the moment the source decays, and the different mass eigenstates propagate with different velocities in the lab frame. Importantly, this is said to lead to an inescapable source of damping of neutrino oscillations at large values of L, since the neutrino wave packets associated with different mass eigenstates separate between production and detection due to their different velocities and may no longer coherently interact with the detector. This is different from our conclusion, that the neutrinos are emitted at different times, such that they arrive at the detector simultaneously, as in Fig. 1. To illustrate the differences between our results and those in the literature: if the source is stationary, maintaining a state unchanged in time, and the particles produced along with the neutrino when the source decays are unmeasured, the common expectation is that oscillations would not persist in the limit L → ∞. However, we find no damping of the oscillations occur in this scenario. In the case of Γ S = 0, the conventional wavepacket picture says that the lighter neutrino wavepacket arrives first at the detector, followed by the heavier neutrino wavepacket. Our calculation shows that the lighter neutrino amplitude is suppressed by e −Γ S ∆t/2 , where ∆t is the difference in transit times, because it is emitted later. Thus the detector sees a "brighter" heavy neutrino emission, simultaneously with a "fainter" light neutrino.
One might think that the oscillation problem has a symmetry given by reversing the direction of time and swapping source and detector. However, in Fig. 1, one fixes the initial state of the quantum system, and evolves it in time. In the time-reversed process, this means the final state rather than the initial state is fixed, which is not the case in experiments. This is an intrinsic asymmetry in the boundary conditions that breaks the time reversal symmetry.
A more realistic treatment of the physics underlying neutrino oscillation experiments can be given using the framework presented here. For example, the neutrino is a fermion, the source has a wave function (though much smaller than any oscillation length), neutrinos are produced and absorbed via the weak interactions, etc. Effects such as these lead to different production and absorption probabilities, but they retain the same oscillation behavior. There are additional physical effects that are present in reactor and pion-beam experiments that can lead to damping of oscillations. These include any treatment of the source and detector that differs from free evolution. In the calculations presented in Sections III and IV, we included the lifetime of the source as a means by which the source does not undergo pure free evolution, but gets exponentially damped in time. Other effects can include small movements of the nucleus in reactor experiments and collisions between the pion and air molecules in accelerator experiments, both of which would lead to damping of the neutrino oscillations. In optics, these are called broadening mechanisms, and have been studied at length (see, for example, Ref. [37,38]). These effects are expected to be negligible for current terrestrial neutrino oscillation experiments. Development of experimental setups where such damping of oscillations is observable would be a valuable verification of the understanding of the quantum mechanical behavior of neutrino oscillations and conceivably could be used to measure previously-unmeasured properties of neutrinos, such as the individual masses of the neutrinos.
The following is in the spirit of the calculation of the behavior of neutrino oscillation presented in the PDG [3]. Consider a single neutrino flavor eigenstate |ν , which is a superposition of mass eigenstates |ν a : On their journey through spacetime, the mass eigenstates all propagate a time T and distance L: and the amplitude for the neutrino to transition from flavor to over time T and distance L is: The oscillation probability is P(T, L) = |A(T, L)| 2 . Because the neutrinos are ultra-relativistic, and if, for simplicity, the neutrino is produced in a 1 → 2 process, S → ν + X, the following approximations can be made for the momentum p a and energy E a for each mass eigenstate: The usual derivation of the neutrino oscillation expression uses the neutrino phase where all neutrino species travel the same distance and the same time T = L, and we have expanded to quadratic order in m a . This equation is the one utilized by experiments when extracting measured values of the neutrino oscillation parameters. The derivation also holds for 1 → many decays, where X is a multiparticle state with invariant mass m X . Note that one does not have to assume "equal energy" nor "equal momentum" among neutrino mass eigenstates, which is a parlance sometimes found in the literature. There are, however, two issues with this standard derivation. First, the justification why one must choose that T = L is not a result of a calculation. One could instead use for T a the neutrino time of flight in the lab frame: This correction in the time of flight can be very small, but including the corrections in Eqs. (A5) and (A6) lead to an oscillation phase that differs from Eq. (A7) in oscillation frequency by a factor of 2: It is not clear which oscillation phase is correct, Eq. (A7) or Eq. (A9), since the only choice was to what order a Taylor expansion is carried out. The second issue is that the oscillation probability P(T, L) depends on the neutrino transit time T . This time is not a measured observable, since there is no information extracted from standard neutrino oscillation experiments regarding when the source particle decays. Therefore, one cannot assign a value to T . Instead, one must sum over the amplitudes for all neutrino transit times. If one does this, as explicitly calculated using a simple system in Sections III and IV, and using a more general system in Appendix B 3, the oscillation phase is unambiguously calculated to be the one in Eq. (A7), and the oscillation amplitude is dominated by terms associated with the source decaying at different times, such that neutrinos of different mass interfere at the detector at the same time, as illustrated with the cartoon in Fig. 1. In neutrino oscillation experiments, the source particle that produces the neutrino must be relatively well-localized in both position and momentum [19]. It is important that the spread in momentum is not zero, so that the source may be localized to a distance much less than the oscillation length: L osc ∆x. As an example, for pion beam experiments with energy of order a GeV: which is much larger than the typical size of the source. The initial state of the position degree of freedom of the source can be approximated as a Gaussian wave function with velocity v centered at x 0 : where d is the number of spatial dimensions, γ is the Lorentz factor γ ≡ (1 − v 2 ) −1/2 , and the subscripts NR indicate non-relativistic state normalization. 6 This state is localized at x 0 with velocity v, and we can make the uncertainties ∆x and ∆v arbitrarily small by taking ∆x → 0 and M → ∞ with M ∆x → ∞. Under free time evolution, the state evolves to 6 The single particle states here have nonrelativistic normalization: where now We can expand the argument of the exponent about p = γM v + k: where we neglect higher-order terms in k, since their effects are suppressed by higher powers of 1/(M ∆x). Performing the Gaussian integral over k, In the regime where t/(M ∆x 2 ) 1, the center of the wave function moves with velocity v and the wave function does not spread with time. Neglecting this term is justified by the fact that e −∆x 2 k 2 suppresses the large-k integration region, so that k ∼ 1/∆x, and the term in question scales like This is negligible if one takes M ∆x 2 → ∞. For experimental sources of neutrinos, e.g., charged pions in accelerator experiments and radioactive nuclei in reactor experiments, the momentum uncertainty ∆x −1 is very small compared to the source's Compton wavelength, as discussed in Section II, and we can ignore the free-particle diffusion of the state |ψ(0) . If so, the wave function simplifies to which, up to a phase, is the same wave function at t = 0, except displaced by an amount vt. It will be useful to note that −iγM (t − v · (x − x 0 )) = −iM t/γ + iγM v · (x − x 0 − vt). If we simultaneously take ∆x → 0, M ∆x → ∞, and M ∆x 2 → ∞, then |ψ(t, x)| 2 can be approximated as a δ-function:

Heavy-to-heavy transitions
Consider the region of parameter space where the source decays to a neutrino and another heavy state, at zero recoil. In this case, we may label the heavy states by an internal label α = G, E, indicating the level of excitation. The energy difference between these two states is the quantity ∆ S appearing in Sections III and IV, i.e., the energy imparted to the daughter neutrino. The regime under discussion here is ∆ S M S , where M S is the mass of the E state of the source. 7 Standard arguments show how, in this regime, the position degree of freedom of the source particle, initially in a state of the form in Eq. (B2), can be treated as a fixed straight-line trajectory, up to corrections 1/M . We review this logic briefly, following Ref. [30]. Consider a scalar theory with fields φ and Φ α , where α = G, E. The action is where X αβ = 0 1 1 0 , and g is a coupling constant. For the moment, we treat φ(x) as a fixed background field that acts as a spacetime-dependent mass for the dynamical field Φ. The two-point Green's function, G αβ (x, y) ≡ Ω| T Φ α (x)Φ β (y) |Ω , satisfies the Schwinger-Dyson equation: The solution to this equation has the path integral representation In a non-relativistic situation, where the spacetime points x and y are separated by a timelike distance large compared to 1/M , the integral is dominated by the stationary-phase configuration, which is a straight-line path from x to y. The exponent in e −iS 1 [x] evaluated on the stationary path x (treating g as a perturbation, and performing the integral over T by the stationary phase as well) is In a sector where there is a single Φ particle, we can simply add this term to the action for the φ field. The first term is just the phase accumulated by the rest energy of the Φ particle. The second term can be rewritten as which is the effective interaction we have used in Sections III and IV.

Heavy-particle states including recoil
If, as a result of the interaction that produces the neutrino, the source particle annihilates into a set of lighter particles, we can no longer treat its trajectory as a straight line throughout the process. Rather, it should be a straight line which ends at the interaction event. This is similar to the situation in HQET of b → u decays, mediated by the current u Γ b v [31]. b v annhilates the b quark moving with velocity v, and u creates the light quarks in the final hadron. Here, we show that the essential form of the amplitude (in particular the phase which leads to the neutrino oscillations) is unchanged as a result of this generalization.
We infer the interaction vertex from a local, translation-invariant interaction vertex between the neutrino 8 field φ(x), and an operator O S→X that annihilates the source state and creates the X state, which is the set of particles produced in conjunction with the neutrino: We have kept the detector infinitely heavy, as in the main text, and used the rotating wave approximation for φ, which is valid for calculations within the light-cone for long baselines. For example, for neutrino production by π → µν µ decay due to the weak Hamiltonian given by projecting onto all π initial states and all µ final states. The S → X transition is given by the form factor F S defined as so that The spin labels on X are included in p X . F S has an additional spinor label α for a fermionic neutrino. Exponential decay of the source can be introduced by the substitution E p S → E p S − iΓ S /2, where 1/Γ S is the lifetime. The leading-order amplitude for neutrino exchange between the source (in the initial state |ψ(t = 0) with momentum space wave function ψ(p S )) and the detector, with an X state of total momentum p X in the final state, is then Using the representation The initial state of the source ψ is given in Section B 1, Eq. (B2). Doing the p S integral, Here, v µ = (γ, γv). This evaluation followed the same logic as in Section B 1. We assumed that F S is analytic in p S , and varies on a typical hadronic scale Λ QCD , which is much bigger than 1/∆x. The source has momentum M v µ , with a small spread of order 1/∆x. Since the form factor is almost constant over this momentum spread, F S can be factored out of the integral Eq. (B23). The remaining integral gives the time-evolved source, which we have seen in Section B 1 moves with velocity v, so that it is now centered at x 0 + vt 2 , as in Eq. (B8). This gives Eq. (B23). ψ(t 2 , x) is the time-evolved spatial wave function: The amplitude is then This is essentially the result derived in the text for a Wilson line source, multiplied by the form factor F S . The significant technical difference is the presence of the wavefunction ψ(x). In the following we show that the integrals can still be done, and the conclusions are unchanged. The impatient reader is invited to skip to the final answer, in Eq. (B61). From now on, we set x 0 = 0 in ψ(t 2 , x), so that the source starts at the origin. Isolating the x dependence of the integral Eq. (B25): where δ ∆p (p) ≡ 1 √ π∆p e −p 2 /∆p 2 is an approximate delta function with width ∆p. Eq. (B25) becomes The exponential ∆ S /γ in Eq. (30) has been replaced by the M/γ − (E X − p X · v), the Doppler shifted energy of the source transition, and the integrand has the form factor F S (p X , M v µ ) and the approximate momentum conserving δ-function δ d 1/∆x (p + p X − γM v). For simplicity, we now set v = 0, x 0 = 0, and Γ S = 0, and proceed to study Doing the t 2 integral, and letting L = Lẑ: Aside #1. It will be useful to estimate the value of the following integral, to leading order in 1/a, 1/b, 1/c, and where a b, c: The integrand is highly peaked at x = 1 if c > 0. Expanding the argument of the exponent around x = 1, where ρ 2 ≡ 1 − x: If instead c < 0, then the integrand is highly peaked around x = −1. Expanding the argument of the exponent around x = −1, where ρ 2 ≡ −1 + x, Using the approximations in Aside #1, one can do the angular integrals, in the limit that ∆x/L 1: The point is that the angular integral over the direction of p, in the limit L ∆x 2 p, is dominated by the saddle point at p =Lp + O ∆x 2 p L -the neutrino comes from the source, up to corrections from the uncertainty in the position of the source.
Changing integration variables from p → E p : Next, deform the contour so that the integration contour is infinitesimally above the real axis: Aside # 2. Consider the following integral: where ∆x, p0 and a are all positive real numbers, and ε is positive and infinitesimally small. The integrand has a pole at E = p0 − iε, and a saddle point at E = p0 − ia/(2∆x 2 ). The goal is to estimate the value of this integral in the limit that a → ∞ and ∆x → ∞. Assuming analyticity, the argument of the exponent can be expressed as the following, letting E = u + iv: Deform the contour away from the positive real axis, going around the pole. Choose a contour C such that the segments are along paths of constant ψ(u, v), and one of the path segments is through the saddle point, as shown in Fig. 3. Then, we have There are two contributions to the integral along C: the location of the saddle point and the location of the integration endpoint. The contribution from the saddle point is If, on the other hand, a < 0, then one can use the same procedure as above, where again the contour has a contribution from the saddle point in the upper half plane and the integration endpoint, but there is no contribution from the pole, since the original contour can be analytically deformed to the upper-half plane: ∆x a e −a 2 /∆x 2 . (B51) The result from Aside #2 implies that the pole is the dominant contribution to the integral, in the limit that ∆x/L 1, which in turn implies that p z X < 0, when t 1 − L is positive (the detector is on completely inside the source's lightcone) yielding: Doing the t 1 integral: A(p X ) −ig S g D (2 3/2 √ π∆x) 3/2 F S (p X , M v µ ) · 1 2(2π) 3 · e i(τ +T /2)(∆ D −M +E X ) 2 sin ((∆ D − M + E X )T /2) Squaring the amplitude, summing over the detector system, and taking there to be only one particle in the X state: Note that the final-state phase space measure lacks the factor of 1/2E X , since we used non-relativistic normalization when we defined in source's wavepacket in Appendix B 1. Using sin 2 (xT /2)/x 2 T 2 πδ(x)+ O(1/T ), do the ∆ D integral, pulling out the form factor F S , since it is the approximately the same for all neutrino species: Let E X = p 2 X + M 2 X . In the limit that m = 0, then the ∆x Gaussian is centered at p * X = − ẑ. In the limit that the support of the ∆x Gaussian is much smaller than the scale of variations of the rest of the integrand, we can expand the argument of the Gaussian around p * X = − ẑ, as well as expanding in m: If the neutrino is ultra-relativistic, the terms proportional to m 2 are very subdominant. So, to a very good approximation, all neutrino species have the same saddle point. Doing the momentum integrals with saddle point methods (here again we assume that the form factor is slowly varying, Λ −1 QCD ∆x, L): × ∞ −∞ dp x dp y dp z e −8∆x 2 M 4 (M 2 +M 2 X ) 2 (pz−p * X ·ẑ) 2 −2∆x 2 (p 2 x +p 2 y ) where E * X ≡ |p * X | 2 + M 2 X . Inputting the expression for C, we have: The conditions for this conclusion to be valid are that L ∆x Λ −1 QCD , using Λ QCD as the typical scale