Consistent QFT description of non-standard neutrino interactions

Neutrino oscillations are precision probes of new physics beyond the Standard Model. Apart from neutrino masses and mixings, they are also sensitive to possible deviations of low-energy interactions between quarks and leptons from the Standard Model predictions. In this paper we develop a systematic description of such non-standard interactions (NSI) in oscillation experiments within the quantum field theory framework. We calculate the event rate and oscillation probability in the presence of general NSI, starting from the effective field theory (EFT) in which new physics modifies the flavor or Lorentz structure of charged-current interactions between leptons and quarks. We also provide the matching between the EFT Wilson coefficients and the widely used simplified quantum-mechanical approach, where new physics is encoded in a set of production and detection NSI parameters. Finally, we discuss the consistency conditions for the standard NSI approach to correctly reproduce the quantum field theory result.

Neutrino oscillations are precision probes of new physics beyond the Standard Model. Apart from neutrino masses and mixings, they are also sensitive to possible deviations of low-energy interactions between quarks and leptons from the Standard Model predictions. In this paper we develop a systematic description of such non-standard interactions (NSI) in oscillation experiments within the quantum field theory framework. We calculate the event rate and oscillation probability in the presence of general NSI, starting from the effective field theory (EFT) in which new physics modifies the flavor or Lorentz structure of charged-current interactions between leptons and quarks. We also provide the matching between the EFT Wilson coefficients and the widely used simplified quantum-mechanical approach, where new physics is encoded in a set of production and detection NSI parameters. Finally, we discuss the consistency conditions for the standard NSI approach to correctly reproduce the quantum field theory result.
Introduction. Precision measurements at low energies are sensitive probes of fundamental interactions, which complement collider searches. Neutrino oscillation experiments [1] are a specific class thereof where one observes a characteristic oscillatory dependence of the neutrino detection rate as a function of the neutrino energy and the distance between the neutrino source and detector. The large body of oscillation data so far has established the existence of at least three distinct neutrino states with different masses [2,3], which is consistent with the predictions of the Standard Model (SM) supplemented with dimension-5 terms leading to Majorana masses for the SM neutrinos [4]. Within this paradigm, the neutrino mass squared differences and the angles of the PMNS mixing matrix have been measured with good accuracy. This opens the door to also probe and constrain new physics (NP), by which we mean NSI between neutrinos and matter that arise from physics beyond the SM (BSM) [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. To this end, however, one needs a map between fundamental parameters of BSM models and observables in oscillation experiments. In this paper we construct such a map for the EFT of the SM degrees of freedom, in which NP modifies the charged-current interactions between neutrinos, charged leptons, and quarks. We also discuss the consistency conditions for the widely used simplified approach, where NP is parametrized by a set of NSI production and detection parameters, to correctly reproduce the quantum field theory (QFT) result.
QFT description. Oscillation probability can be rigorously derived in the framework of quantum field theory. Various derivations are available in the literature in the absence of NSI (see e.g. [25][26][27]). Below we give an expression valid for completely general interactions between neutrinos and matter. Consider neutrinos produced in the process S → X α ν (e.g. beta decay of a nucleus in a reactor, or pion decay), where X α is one or more body final states containing one of the charged leptons α = (e, µ, τ ). The neutrinos are detected via the process ν T → Y β (e.g. inverse beta decay), where again Y β contains a charged lepton β . The production and detection can be described by QFT amplitudes , where the index k labels neutrino mass eigenstates. The information about fundamental parameters, is encoded in M P αk and M D βk , which should be then connected to observables. For the source (S) and target (T) states separated by a macroscopic distance L, the observable is the differential rate of detected events R αβ ≡ dN αβ /dtdE ν given by . A compact derivation of this formula is presented in Appendix A, where we also enumerate its limitations. Above, ∆m 2 kl ≡ m 2 k − m 2 l is the mass squared difference between the neutrino eigenstates. The phase space elements dΠ P,D for the production and detection processes are defined in the standard way: where P is the total 4-momentum of the initial state and k i are the 4-momenta of the final states. For the production dΠ P includes the neutrino phase space d 3 kν (2π) 3 2Eν and we define dΠ P via dΠ P ≡ dΠ P dE ν . The sign in Eq. (1) involves integration and sum/average over all unobserved degrees of freedom, such as angular variables and spins. Finally, complex conjugated amplitudes are denoted with a bar, N S,T are the number of source and target particles, m S,T are their masses, and E ν is the neutrino energy, which is an observable in typical experiments.
The rate in Eq. (1) displays the famous oscillatory be- factor. In the absence of oscillations, rates would be calculated using the neu-arXiv:1910.02971v2 [hep-ph] 29 Oct 2019 trino differential flux Φ α ≡ N S N T 4πL 2 dΓ P α dtdEν , and the detection cross section at the target. The source decay rate Γ P α (with an emission of α and summed over neutrino mass eigenstates) and the detection cross section σ D β (with an emission of β and summed over neutrino mass eigenstates) can be calculated by the usual means in QFT. We have One can define the ν α → ν β oscillation probability as the ratio of the rate of detected events in Eq. (1) to the no-oscillation expression in Eq. (2), finding This formula appears in Ref. [28] in a slightly different form without explicit phase space integration. Oscillation probability is an intuitive and widely employed concept, however strictly speaking P αβ is not an observable. For this reason in the rest of this paper we work with the rate in Eq. (1). QM-NSI description. The machinery of QFT is rarely employed in the neutrino literature. Instead, a simpler quantum mechanical (QM) description of neutrino oscillations is most often used. One defines flavor states as linear combinations of mass eigenstates: |ν α = k U αk |ν k , where U is the unitary PMNS mixing matrix. In this language, the NSI effects are encoded in parameters s,d αβ , which correspond to non-standard effects in neutrino production and detection, respectively [19,29,30]. They are defined by the mismatch between pure flavor states and the neutrino states produced at the source and detected at the target, namely [9]: with the normalization N s For anti-neutrinos Eq. (4) holds with s,d → ( s,d ) * . The probability of |ν s α oscillating into |ν d β is given by where in the absence of the matter effects in the propagation the Hamiltonian is H βα = k U βk m 2 k U * αk /(2E ν ). In this approach, which we refer to as the QM-NSI formalism, the event rate is then given by [6] where x s ≡ (1 + s )U * and x d ≡ (1 + d ) T U . Above, Φ SM α and σ SM β are the incident flux and detection cross section calculated in the absence of NSI. The normalization factors N s α N d β cancel in the observable rate and thus one could have omitted them altogether [6]; their only role is to ensure that P αβ ≤ 1, that is it can be interpreted as a probability.
Results from oscillation experiments are often presented or recast as constraints on the NSI parameters s,d αβ . However, the utility of the latter hinges on whether they can be unambiguously connected to more fundamental parameters of the microscopic theory. Only after such matching the coefficients s,d αβ determined in different experimental settings can be meaningfully compared and combined. In the following we discuss this issue, and illustrate it with physically relevant examples. We will define the conditions under which the NSI parameters can indeed provide an adequate description of NP effects in neutrino oscillation. Conversely, we will show examples where this is not the case.
Matching QFT and QM-NSI results. One could try to match the QM-NSI and QFT language starting from the definition in Eq. (4). This however would be problematic, as such concepts as neutrino flavor states or production and detection states are murky in a QFT framework when general neutrino interactions with matter are allowed. Therefore, we will follow a pragmatic approach and match the observable rates predicted by the QFT (Eq. (1)) and QM-NSI frameworks (Eq. (5)). This comparison will allow us to determine the map between the NSI s,d and the Lagrangian parameters, or else conclude the map does not exist.
In this paper we focus on NP in charged current interactions between neutrino and matter. The microscopic theory we consider is the EFT of the SM degrees of freedom at the energy scale µ ≈ 2 GeV, in which NP modifies the effective 4-fermion interactions between leptons and quarks (extensions to other theories and interactions are straightforward using our approach). At leading order in this EFT the neutrino interactions with matter can be parametrized by the Lagrangian where and P L,R are the usual chirality projectors (1 ∓ γ 5 )/2. The quarks u, d, and charged leptons α are in the basis where their kinetic and mass terms are diagonal. For the neutrino fields the kinetic terms are diagonal but the mass terms are not: they are connected to the mass eigenstates by the unitary rotation via the PMNS matrix, ν α = k U αk ν k . In this EFT the effects of NP are parametrized by the Wilson coefficients [ X ] αβ , which encode new interactions between quarks and leptons mediated by BSM particles heavier than ∼ 2 GeV. For example, non-zero R can arise in left-right symmetric models due to the W boson coupling to right-handed quarks and mixing with the SM W , while non-zero S,P,T are generally predicted in leptoquarks models. More generally, X can be connected to parameters of the weak-scale EFT, known as the SMEFT [31][32][33].
SM interactions. To warm up, let us first calculate the oscillation probability in the limit of SM interactions, which corresponds to setting all X = 0. In this case, which was studied in Ref. [25], the amplitudes can be decomposed as M are independent of the neutrino mass index k up to negligible corrections, whereas they do depend on the charged-lepton flavor index α (which we omit to ease the notation). They also depend on the kinematic and spin variables in the production and detection processes, and they appear in the observables integrated/averaged over by dΠ P ,D . All in all the rate in Eq. (1) can be written as where the SM flux and cross-section are given by Exactly the same result is obtained from Eq. (5) in the limit s,d αβ = 0. V -A interactions. A less trivial example is when NP enters only via V -A interactions: [ L ] αβ = 0 [5,6,34,35]. In this case the detection/production amplitudes decompose as M In fact, the quantity x L is equivalent to a "non-unitary mixing matrix", an approach that has been studied in neutrino literature [6,35]. The same result is obtained in the QM-NSI approach from Eq. (5) when the NSI parameters are mapped to the Lagrangian parameters as [35] V−A : In the V -A case the map between NSI and Lagrangian parameters is well-defined, unambiguous, and simple. The NSI parameters for production and detection are the same up to Hermitian conjugation. General case. For general NP interactions in Eq. (6), the production and detection amplitudes can be decomposed as The sum above goes over all types of interactions in Eq. (6): X = L, R, S, P, T . We stress that A P,D X will typically have completely different dependence on kinematic and spin variables for different X. Plugging this decomposition into Eq. (1) we obtain a lengthy expression, which we nevertheless quote in full: where we define the production and detection coefficients We show in Appendix B the expressions of the above coefficients for different processes. For anti-neutrinos Eq. (12) holds with U ↔ U * and X ↔ * X . The formulas for the neutrino oscillation probability are collected in Appendix C.
At the linear level in the QFT expression in Eq. (12) matches the QM-NSI one in Eq. (5) provided the NSI parameters are expressed by the EFT parameters as Therefore the QM-NSI formalism can approximate the correct oscillation probability obtained from the general EFT as long as the deviation from the SM, encoded in the coefficients [ X ] αβ , is sufficiently small. If non-V -A interactions are involved, the NSI parameters obtained via the matching in Eq. (14) may be a function of the neutrino energy. The production and detection parameters depend on the same X parameters, but they do not satisfy anymore the relation s = d † valid for the V -A case.
Beyond the linear approximation the QM-NSI formalism fails in general because no matching can be found to connect with the QFT result. This is one of our main results. The consistency condition for the matching in  Eq. (14) to be valid to all orders in is for all X and Y for which X,Y are non-zero. Eq. (15) is trivially satisfied if the only NP deformations are of the V -A type, that is if only L is non-zero, in agreement with our previous discussion. However, for non-V -A deformations Eq. (15) is typically not satisfied, because then A P,D X may have different dependence on kinematic variables than A P,D L . We now give a concrete example where Eq. (15) is not satisfied. In many oscillation experiments low-energy neutrinos are detected via inverse beta decay ν p → n e. Consider NP of the V +A type affecting the process, [ R ] eβ = 0 for some β. In this setting we find d RL = , and d RR = 1 where g A is the axial charge of the nucleon. Thus we see that |d RL | 2 = d RR , since g A = 1.251 (33) [36][37][38]. We conclude that the effect of V +A interactions in neutrino experiments that involve inverse beta decay cannot be described by the NSI parameters s,d beyond the linear level. In the presence of scalar and tensor interactions the detection coefficients are given in Eq. (B3), and again |d SL | 2 = d SS , |d T L | 2 = d T T . Moreover, in those two cases the left-hand sides depend on the neutrino energy, while the right-hand sides do not. We show in Fig. 1 an example of a mismatch between the rates calculated in the QM-NSI and QFT formalisms.
One should not however jump into conclusion that the QM-NSI formalism always fails for non-V -A interactions. For instance, this is not the case if A P,D X = c P,D X A P,D L , where c P,D X is a constant independent of the kinematic variables to be integrated over. The latter happens e.g. in the 2-body decay of a spin-zero particle. This includes of course the phenomenologically relevant case of neutrino production through pion decays. Thanks to the pseudoscalar nature of the pion, the only non-zero hadronic matrix elements for this decay are 0|ūγ µ γ 5 d |π + , and 0|ūγ 5 |π + . As a result the production process is sensitive only to axial ( L -R ) and pseudo-scalar ( P ) interactions. We list all the non-zero production coefficients in Eq. (B5). We see that for all non-V -A interactions the consistency condition is satisfied: p XL p * Y L = p XY for X, Y = R, P . Therefore, neutrino production via pion decays can be described by the QM-NSI formalism to all orders, even in the presence of non-V -A interactions. In Table I we summarize the linear matching between the NSI parameters and EFT Wilson coefficients for several relevant production and detection processes.
Discussion and conclusions. We close with several comments on the results derived above: 1. In this paper we only discussed charged-current NSI and assumed the absence of matter effects in propagation. The neutral-current NSI other than the matter effects can also be correctly described by QFT expressions analogous to Eqs. (1)-(3), and they are relevant e.g. if neutrinos are detected via coherent scattering on nuclei. To include NSI entering via the matter effects one would need to modify the neutrino propagator in the derivation in Appendix A starting from Eq. (A6). We leave this for future work.
2. It is worth stressing that charged-current NSI modify not only the flux and cross-section in Eq. (2), but also the oscillation probability in vacuum. The latter follows directly from Eq. (3), which depends on the production and detection amplitudes. Generically, that dependence does not cancel between the numerator and denominator in Eq. (3).
3. The observable in Eq. (1) may depend on NP in two distinct ways. One is direct, e.g. through a dependence of the production and detection amplitudes M P,D αk on the NP parameters X of the Lagrangian in Eq. (6). The other is indirect, due to NP "polluting" the observable used for determination of the SM parameters [39]. This is the case for the CKM parameter V ud in Eq. (6). If NP is present, β decay experiments determine a combination of V ud and [ X ] eβ parameters, and in this case the value of V ud cannot be just taken from PDG. This indirect effect is ignored in most of the prior neutrino literature, even though it is of the same order as the direct effects, leading to incorrect results. For instance, indirect and direct effects generated by the coefficient [ X ] ee cancel at all orders, making this coefficient unobservable in oscillation experiments [32].   lation probability that go beyond its obvious intuitive qualities. First, for the sake of calculating ratios of measurements at different distances L for a fixed E ν , the ratio of probabilities is the same as the ratio of rates. In fact, many neutrino analyses use only the ratios to bypass large uncertainties in overall neutrino rates, and for this purpose the oscillation probability is sufficient. Second, in many familiar scenarios the physics contributing to the oscillation probability in Eq. (3) and no-oscillation piece in Eq. (2) is distinct. This is the case in the SM, with electroweak and hadronic parameters contributing to the flux and cross-section, and neutrino masses and mixing angles contributing to the oscillation probability. Another kind of separation happens in the EFT scenario (6) at the linear level, where flavor off-diagonal X affect the oscillation probability, whereas flavor diagonal X affect the flux/cross section [32]. It is important to note however that such a separation does not hold in a general BSM scenario. In particular, quadratic effects of both diagonal and off-diagonal X contribute to both pieces.

5.
We have not really defined the basis for the neutrino states in the effective Lagrangian in Eq. (6). A natural flavor basis exists if that Lagrangian is derived from the underlying SMEFT theory with manifest electroweak symmetry. Namely, starting from a generic basis with lepton doublets L α = ( α , ν α ), one rotates the doublets to the basis where the charged leptons mass matrix is diagonalized, L α = ( α , ν α ) = V e L α . In this language, the PMNS matrix describes the mismatch between charged and neutral lepton rotations to the mass basis, in analogy to the CKM matrix in the quark sector. However, from the viewpoint of the EFT below the weak scale, where charged and neutral leptons are not collected in doublets, nothing distinguishes the basis of ν α in Eq. (6) as soon as [ L ] αβ = 0. Unitary rotations ν α → V ν α transform the Lagrangian into another equivalent form with NP parameters rotated as X → X V and the neutrino mass matrix rotated by M ν → V T M ν V . Physics of course cannot depend on which of the infinite family of bases we choose to work with. And indeed, the observable rate in Eq. (12) is invariant under unitarity rotations X → X V accompanied by U → V † U .
We find that there are no such "zero-distance effects" at linear order. Let us note that in the α = β case the rate itself is affected by linear effects in [ X ] αα , but they come from NP modifying the neutrino flux and detection cross-section in Eq. (2). At quadratic order, zero-distance effects do appear in general. On the other hand, they vanish at all orders in the α = β case with V -A interactions, i.e. P V −A αα (L = 0) = 1. Our results are relevant for the study of zero-distance effects since they are quadratic and, in the α = β case, necessarily non-V -A.
In conclusion, the main results of this paper are: i) The expression in Eq. (12) for the event rate in neutrino oscillation experiments including nonstandard chargedcurrent interactions described by the EFT Lagrangian in Eq. (6); ii) The matching in Eq. (14), valid at the linear level in NP, between the EFT coefficients that describe the underlying interactions and the QM-NSI parameters; iii) The consistency condition in Eq. (15) for that matching (and the simplified QM-NSI approach itself) to be valid to all orders in NP parameters.
Our results are particularly relevant for analyses of oscillation data when effects of non-V -A physics (or equivalently s = d † ) are taken into account [7,9,[11][12][13][15][16][17]. We clarify the microscopic meaning and validity of the long list of existing analyses of oscillation data carried out within the traditional QM-NSI approach. Their discovery potential can now be consistently analyzed and compared, among themselves and together with nonoscillation probes that are sensitive to the same underlying physics. In this appendix we derive the master formula in Eq. (1) describing the number of neutrino events detected at a distance L from the source, taking into account possible neutrino oscillations and nonstandard charged-current interactions. Our approach follows similar steps as Ref. [25]. The two main differences are: 1) we allow for general non-SM charged-current interactions in neutrino production and detection, and 2) we work with time-independent packets for the source and target particles, which greatly simplifies further mathematical transformations. Of course, the source is necessarily unstable, hence the latter assumption will lead to one subtlety in the derivation below.
We consider an experimental setup where neutrinos are produced in a process A x → X α ν, and detected via the process νB y → Y β . Here X α and Y β are n x -and n y -body final states (n i ≥ 1). The indices α and β indicate that these states contain charged lepton α and β respectively, but otherwise their precise identity is irrelevant for this discussion. A x and B y are both one-body particle states localized in the coordinate space, describing the neutrino source (e.g. a beta-decaying nucleus in a reactor) and target (e.g. a proton in a detector). We will work in the time-independent approximation where the states A x and B y are represented by wave-packets of scattering in-states which do not change in time. We parametrize them as where E j = m 2 j + |p j | 2 , for j = A, B, | p j in are momentum eigenstates normalized as q j | p j = (2π) 3 2E j δ 3 ( p j − q j ), and the states are normalized as A x |A x = B y |B y = 1. For simplicity we choose Gaussian wave packets for both states with the same spread σ in the position space: The wave packet describes a particle at rest localized near z with the uncertainty of order σ.
The idea is to treat the neutrino production and detection together as a single process [25]: rather than consider the neutrino production and detection separately. In this approach, neutrino is merely an intermediate particle in the amplitude. The outgoing states are approximated by pure momentum eigenstates with the eigenvalues k i , i = 1 . . . n, where n = n x + n y is the number of particles in the final states. We are interested in the transition probability for this process: Plugging the wave packets for the initial states, and using out k 1 k 2 . . . k n |p A p B in = (2π) 4 δ 4 (p A + p B − k i )M we obtain Up to this point, we followed the classic derivation of the cross section formula, see e.g. Ref. [41]. What distinguishes the case at hand is the particular choice of the initial states A x |, B y | describing two spatially separated particles (rather than head-on beams as in the cross section case). Furthermore, the amplitude for this process is dominated by the kinematic region where q ≡ p A − p X = p Y − p B is close to the neutrino mass shell, q 2 ≈ m 2 k . In that region, unitarity requires the amplitude to factorize into the production, detection, and propagation parts: where the sum goes over all neutrino eigenstates, and the amplitudes in the numerator are evaluated for all particles on-shell, including the neutrinos. We can also factorize the phase space: dΠ n = dq 2 0 2π dΠ P dΠ D , where the first factor is the X+neutrino phase space, and the second factor is the Y phase space. Next, it is convenient to isolate the phases in the wave packets by rewriting φ z (p) =φ(p)e i p z . Finally, we trade one delta function for a time integral where L = y − x. The next step is to perform the q 2 0 integral, treating it as a contour integral: where ∆m 2 kl = m 2 k − m 2 l . Above, the amplitudes are now evaluated at q 0 = | q| 2 + m 2 k , that is for on-shell neutrinos. At this point we introduce a number of approximations that are appropriate for the description of broad classes of real-life neutrino experiments: 1. The intermediate neutrinos are relativistic, hence in the production and detection amplitudes we can set q 0 = | q|.
The dependence on the neutrino masses survives only via the ∆m 2 kl factor in Eq. (A8). 2. The wave packets are localized in an area much larger than the inverse mass of the source and target particles, In the subsequent analysis we will only keep terms of O(σ −1 ) and ignore those of O(σ −2 ). In particular, we can approximate E A ≈ m A and E B ≈ m B .
3. We ignore the dependence of the amplitudes on p j or p j , from which it follows that M = M. Given the assumption in the previous point, this present assumption is safe whenever the amplitudes are dominated by a velocity-independent term.
With the above assumptions Eq. (A8) simplifies to Due to our approximations, after replacing e i(E A +E B −E A −E B )t ≈ 1 nothing in the integrand depends on t and thus the integral is infinite. This singularity could in fact be expected: due to using time-independent wave packets for the source |A x we tacitly integrate the rate of the A x B y → XY process from t = −∞ to t = +∞. In a physical situation, however, A x is unstable, appears at a finite time t 0 , and decays after a finite time t 0 + T . Outside this window the process A x B y → XY cannot occur. With this in mind, we drop the integration over t, and obtain the following result for the rate, that is the number of events per unit time: Next, it is convenient to change the integration variables as p ± j = p j ± p j . Afterwards we can trivially perform the Gaussian integral over d 3 p + A d 3 p + B ,and eliminate the integral over d 3 p − A using the δ 3 . We also fix the coordinate frame such that L = (0, 0, L), so that the z-axis connects the source and the target. This leads to where we simplified the notation: p ≡ p − B = − p − A . In the integration over p z , the principal value is suppressed by the rapidly oscillating e iLpz , and is neglected, which leaves the contribution from the pole at p z = (∆m 2 kl − 2q i p i )/2q z : where i = x, y. Note that q i is the neutrino momentum in the "wrong", off-axis direction transverse to L, thus |q i | |q z | as long as L σ. Therefore we can drop q i /q z factors everywhere, except when they are multiplied by L. Then we can trivially perform the Gaussian integral over p x and p y , which leads to (A13) The oscillatory factor e −i2πL/Losc appears here for the first time in this derivation, with the oscillation length L osc = 4πqz ∆m 2 kl ≈ 4πEν ∆m 2 kl . In the QFT approach it arises because of interference between distinct neutrino mass eigenstates k = l entering the propagators in Eq. (A7), which in turn is possible due to momentum uncertainty described by the initial state wave packets. The oscillations become suppressed by the factor e −∆m 4 kl σ 2 /8q 2 z when the packet size becomes comparable to the oscillation length [42]. A condition for oscillations to be possible is On the other hand, in our approach we do not find exponential suppression of the oscillations proportional to the distance L travelled by the neutrino. The usual argument for this suppression [42,43], due the decoherence of wave packets corresponding to different neutrino eigenstates traveling at different speeds, does not apply in the static situation considered here. Due to the last exponential factor in Eq. (A13) the neutrino production angle θ ≈ q 2 x + q 2 y /q z must be such that Lθ σ. This has a simple physical interpretation: the neutrino must hit the target within its position uncertainty described by the wave packet. Neutrinos emitted with Lθ σ simply miss the target and do not contribute to the probability of the A x B y → XY process. With this in mind, on the final transformation we trade q z = E ν cos θ ≈ E ν , and q 2 x + q 2 y = E 2 ν sin 2 θ ≈ E 2 ν θ 2 . The production phase space contains the neutrino phase space d 3 q (2π) 3 2q0 = Eν dEν d cos θdφ 16π 3 . One more assumption we need is that neutrinos are produced isotropically, that is M P αkM P αl integrate/summed over unobserved degrees of freedom is independent of the angular variables θ, φ. This assumption is satisfied in typical neutrino experiments where the source is unpolarized. The integral over θ can be evaluated order by order in σ 2 /L 2 , leading to where dΠ P = dΠ P dE ν . Note the geometric 1/L 2 factor in front, which is of course expected intuitively. Mathematically, it appears due to integrating over the neutrino production angles in the phase space, where the contribution of off-axis neutrinos is exponentially suppressed and only the small cone θ σ/L effectively contributes to the transition rate. The dependence on the size σ of the initial wave packets has canceled out, except in the last exponential factor, which can be ignored in the limit σ∆m 2 kl E ν . In that limit, after multiplying the rate by the number of source and detector particles N S,T we obtain the master formula in Eq. (1).