text:
# Towards a unique formula for neutrino oscillations in vacuum

## 1 Introduction

Nowadays the study of neutrino oscillations is the particle physicists’ best hope to learn how to extend the Standard Model, since it provides information on the neutrino masses and mixings which are thought to originate in new physics. Neutrino oscillations are an example of flavor oscillations, which arise when the particles produced and detected in an experiment are superpositions of different mass eigenstates. The interference between the propagating mass eigenstates leads to a spatial oscillation of the detection probability of a neutrino with a given flavor. The standard oscillation formula, expressing the probability to detect a neutrino of momentum $`p_0`$ emitted with a flavor $`\alpha `$ as having a flavor $`\beta `$ after a propagation on a distance $`L`$, reads

$$𝒫_{\alpha \beta }(L)=\underset{i,j}{}V_{i\alpha }V_{i\beta }^{}V_{j\alpha }^{}V_{j\beta }\mathrm{exp}\left(i\frac{m_i^2m_j^2}{2p_0}L\right),$$
(1)
where $`m_{i,j}`$ are the neutrino mass eigenstates and $`V`$ is the matrix relating the flavor to the mass fields. This phenomenon can explain the depletion of the solar electron-neutrino flux , as well as the up/down asymmetry of the neutrino flux produced by cosmic rays in the atmosphere . Moreover, the long baseline experiment K2K hints to a muon-neutrino depletion on a terrestrial distance , while controversial evidence for neutrino oscillations has appeared in the laboratory experiment LSND .

In view of the importance of neutrino oscillations, it is rather unsettling that the neutrino oscillation formula (1) used to analyze experimental results is nearly always derived in a completely inconsistent way. In particular, this derivation usually starts with unlocalized plane wave states and ends, strangely, with a space-dependent oscillation formula. This hocus-pocus requires at least two arbitrary assumptions: the equal time prescription, stating that the propagation time is the same for the different mass eigenstates, and the classical propagation condition, $`x=vt`$ (see  for a review). Additional assumptions have generated endless discussions about the equality of the energy-momenta of the different mass eigenstates. Moreover, the rejection of the equal time prescription by some authors has lead them to predict oscillation lengths differing from the standard result by some factor (a factor 2 is often quoted).

A more consistent derivation of the neutrino oscillation formula has been done in the intermediate wave packet model, in which the neutrino mass eigenstates are represented by wave packets propagating in space-time . Although the classical propagation condition can be dispensed with, this method implicitly requires the use of the arbitrary equal time prescription in order to obtain the standard result for the oscillation length. Moreover, this derivation is not wholly consistent because it involves flavor eigenstates which cannot be well-defined . Finally, while the 3-momentum uncertainty is naturally included in this model, the energy uncertainty has to be put in by hand, with the result that there is no agreement on whether the coherence length has an upper bound or not .

In order to solve the various problems of the intermediate wave packet model, different models using quantum field theory have been proposed. They can be grouped in four categories. The first category groups the external wave packet models, the best example of which is given by the seminal paper by Giunti, Kim, Lee and Lee . In these models, the neutrino is considered as an internal line of a Feynman diagram, propagating between a source and a detector, which are represented by in- and outgoing wave packets . A variation on this theme is the Kiers-Weiss model, in which the external wave packets are replaced by quantum oscillators . The second category groups the stationary boundary conditions models, the best example of which is given by the Grimus-Stockinger model . They are very similar to the external wave packet models, except that the wave packets are replaced by states with a well-defined energy, i.e. stationary states . The third category groups models where the neutrino is represented by its propagator and coupled at production with a source, though not coupled at detection with anything . Note that this kind of model requires the equal time prescription, as in the intermediate wave packet model, in order to obtain the standard oscillation length, since the detection mechanism is left unspecified. The fourth category includes the Blasone-Vitiello model, in which the construction of a Fock space of flavor eigenstates is attempted .

While few authors deny that the most rigorous treatment of oscillations is done in a quantum field framework, the variety of such models has not made a good case for their widespread acceptance. Besides, these models sometimes yield conflicting oscillation formulas, either because of their different assumptions or because of their different approximation methods. The main aim of this article is to show which of these results are correct and how they can all be understood within the external wave packet model. In particular, the Grimus-Stockinger oscillation formula will be seen to be a limiting case of the Giunti-Kim-Lee formula, while the Ioannisian-Pilaftsis and Shtanov new oscillation formulas will be shown to be false. At the same time, the role of accurate energy measurements in the coherence length will be clarified.

Contrary to the existing literature, our calculations take into account the dispersion and do not resort to a relativistic limit. In this way, the final oscillation formula will be valid as well for neutrinos as for K and B mesons (it will be shown that the spin can be neglected). This unification of the neutrino and meson oscillation formulas makes it much more difficult to accept nonstandard results, since the meson oscillations are well studied and the meson oscillation length has been cross-checked with different methods. Another of our concerns will be to underline the quantum mechanical interpretation of the external wave packet model. The quantum field theory models are often perceived as complicated and obscure, whereas they also exist in simple versions  and have a clear physical interpretation in direct correspondence with intermediate wave packet models.

The outline of the article is as follows. In the second section, the gaussian external wave packet model is reviewed and all useful definitions are given. The flavor-mixing amplitude is expressed as a convolution of the neutrino propagator with a function depending on the overlap of the external wave packets. This section does not contain new material and follows the notation of . In the third section, a new way to evaluate the convolution is explained, leading to an explicit oscillation formula. It is then shown that the coherence length, beyond which oscillations vanish, can be increased without bound by more accurate energy measurements. As a corollary, it is proved that the stationary boundary conditions model is a special case of the external wave packet model. In the fourth section, another new method to evaluate the convolution is presented, yielding the explicit dependence of the amplitude on time and distance. This method allows a clear physical interpretation of the amplitude in terms of wave packets. In particular, the propagation range is subdivided into three regimes by two thresholds marking the onset of the transversal and longitudinal dispersion of the wave packet. With this method, the new oscillation formulas proposed in Refs.  and  can be disproved. Finally, we show that the Blasone-Vitiello oscillation formula  is phenomenologically equivalent to the standard result.

## 2 The external wave packet model

The propagating process of a particle between a source and a detector (indicated by dotted circles) is symbolized by the following diagram:

FIG. 1. Propagation of a particle $`\nu `$ between a source and a detector, centered in $`x_{\mathrm{SS}P}`$ and $`x_{\mathrm{SS}D}`$.

The arrows indicate the momentum flow. $`P_I`$ represents the set of incoming particles, of total momentum $`q`$, arriving in the production region (or source), which is centered around the point $`(t_{\mathrm{SS}P},𝐱_{\mathrm{SS}P})`$. $`P_F`$ represents the set of outgoing particles, of total momentum $`k`$, coming from the production region, with the exception of the intermediate neutrino $`\nu `$ whose propagation is studied. $`D_I`$, $`D_F`$ and $`(t_{\mathrm{SS}D},𝐱_{\mathrm{SS}D})`$ are defined similarly, but apply to the detection process. The interaction points at production and detection are noted $`x`$ and $`x^{}`$, respectively. All particles are assumed to be stable. The internal line could also represent an antineutrino, but we shall assume that the experimental conditions are such that a quasi-real neutrino propagates on a macroscopic distance, transferring positive energy from $`x`$ to $`x^{}`$. If $`x^0x^0`$ is a macroscopic time, it will be seen that the neutrino $`\nu `$ contributes to the propagation, but not its antiparticle.

Since known neutrinos are much lighter than the experimental energy thresholds, their masses can be neglected in the interaction vertices and in the numerator of the neutrino propagator. In that case, the spin structure factorizes from the sum over the mass eigenstates, so that the computation of the oscillation formula can be carried out as if the neutrinos were scalar. If one wants to consider nonrelativistic neutrinos, factorization also occurs for nearly degenerate masses. The only case where factorization does not occur is nonrelativistic neutrinos with very different masses. However, spatial oscillations vanish in these conditions since it is experimentally possible to determine which mass eigenstate contributes to the amplitude. In that case, a treatment with wave packets is not necessary, since the flavor-changing probability does not depend on the propagation distance. In the following we shall assume that the neutrinos are either relativistic or nearly degenerate in mass. Thus the spin structure can be absorbed in an overall multiplying factor.

The external wave packets are expressed in terms of their momentum-space wave function: $`|\psi >=[d𝐤]\psi (𝐤,𝐊)|𝐤>`$, where $`[d𝐤]`$ is the integration measure in 3-momentum space and $`𝐊`$ is the average 3-momentum. If $`\psi (𝐤,𝐊)`$ is the momentum space wave function of a wave packet centered in $`𝐱=\mathrm{𝟎}`$ at time $`t=0`$, then

$$\mathrm{\Psi }(𝐤,𝐊,𝐱_0,t_0)=\psi (𝐤,𝐊)e^{iE(𝐤)t_0i𝐤𝐱_0},$$

is the momentum-space wave function of a wave packet centered in $`𝐱_0`$ at time $`t_0`$.

Without loss of generality, let us choose to work with only one particle in $`P_I(q)`$, in $`P_F(k)`$, in $`D_I(q^{})`$ and in $`D_F(k^{})`$. The extension to a larger number is straightforward and would only complicate the notation. The wave packets are built such that those involved in the production of the $`\nu `$ are centered at $`𝐱_{\mathrm{SS}P}`$ at time $`t_{\mathrm{SS}P}`$, while those involved in the detection of the $`\nu `$ are centered at $`𝐱_{\mathrm{SS}D}`$ at time $`t_{\mathrm{SS}D}`$. They are noted

$`|P_I>`$ $`=`$ $`{\displaystyle [d𝐪]\mathrm{\Psi }_{P_{\mathrm{SS}I}}(𝐪,𝐐,𝐱_{\mathrm{SS}P},t_{\mathrm{SS}P})|P_{\mathrm{SS}I}(𝐪)>}`$
$`|P_F>`$ $`=`$ $`{\displaystyle [d𝐤]\mathrm{\Psi }_{P_{\mathrm{SS}F}}(𝐤,𝐊,𝐱_{\mathrm{SS}P},t_{\mathrm{SS}P})|P_{\mathrm{SS}F}(𝐤)>}`$
$`|D_I>`$ $`=`$ $`{\displaystyle [d𝐪^{}]\mathrm{\Psi }_{D_{\mathrm{SS}I}}(𝐪^{},𝐐^{},𝐱_{\mathrm{SS}D},t_{\mathrm{SS}D})|D_{\mathrm{SS}I}(𝐪^{})>}`$
$`|D_F>`$ $`=`$ $`{\displaystyle [d𝐤^{}]\mathrm{\Psi }_{D_{\mathrm{SS}F}}(𝐤^{},𝐊^{},𝐱_{\mathrm{SS}D},t_{\mathrm{SS}D})|D_{\mathrm{SS}F}(𝐤^{})>}.`$
Let us first suppose that there is only one mass eigenstate $`m_j`$ associated to the propagating neutrino of FIG. 1. The general formula of the connected amplitude corresponding to this process is given by

$$𝒜_j=<P_F,D_F|\widehat{T}\left(\mathrm{exp}\left(id^4x_I\right)\right)\mathrm{𝟏}|P_I,D_I>,$$

where $`_I`$ is the interaction Lagrangian for the neutrino and $`\widehat{T}`$ is the time ordering operator. Let $`g`$ be the coupling constant of $`\nu `$ with the other fields. Expanding the amplitude to order $`g^2`$, and inserting the wave packets expressions, one obtains

$$𝒜_j=i\frac{d^4p}{(2\pi )^4}\frac{\psi (p^0,𝐩)}{p^2m_j^2+iϵ}e^{ip^0T+i𝐩𝐋},$$
(2)
where the average propagation time $`T`$ is defined by $`T=x_{\mathrm{SS}D}^0x_{\mathrm{SS}P}^0`$ and the average propagation distance by $`𝐋=𝐱_{\mathrm{SS}D}𝐱_{\mathrm{SS}P}`$.

The overlap function $`\psi (p^0,𝐩)`$ represents the overlap of the incoming and outgoing wave packets, both at the source and at the detector. It is defined by

$`\psi (p^0,𝐩)={\displaystyle d^4xe^{ipx}d^4x^{}e^{ipx^{}}[d𝐪]\psi _{\mathrm{SS}P_{in}}(𝐪,𝐐)e^{iqx}[d𝐤]\psi _{\mathrm{SS}P_{out}}^{}(𝐤,𝐊)e^{ikx}}`$
$`\times {\displaystyle }[d𝐪^{}]\psi _{\mathrm{SS}D_{in}}(𝐪^{},𝐐^{})e^{iq^{}x^{}}{\displaystyle }[d𝐤^{}]\psi _{\mathrm{SS}D_{out}}^{}(𝐤^{},𝐊^{})e^{ik^{}x^{}}M_P(q,k)M_D(q^{},k^{}),`$ (3)
where $`M_P(q,k)`$ and $`M_D(q^{},k^{})`$ are the interaction amplitudes at production and detection. All external particles are on-shell, i.e. $`q^0=E_{P_{\mathrm{SS}I}}(𝐪)=\sqrt{𝐪^2+m_{P_{\mathrm{SS}I}}^2},`$ and so on. Note that the overlap function is independent of $`x_{\mathrm{SS}P}`$ and $`x_{\mathrm{SS}D}`$. The integrals over $`x`$ and $`x^{}`$ in Eq. (3) yield delta functions, which impose energy-momentum conservation at the source and the detector. The amplitude corresponding to the propagation of an antineutrino from $`x_{\mathrm{SS}P}`$ to $`x_{\mathrm{SS}D}`$ can be obtained by the substitution $`pp`$.

Let us now consider the case where the intermediate state in FIG. 1 is a superposition of different mass eigenstates. The fields in the flavor and mass bases are related by a unitary matrix $`V`$

$$\nu _\alpha =\underset{j}{}V_{\alpha j}^{}\nu _j,$$

where Greek and Latin indices refer to the flavor and to the mass basis, respectively. We assume that it is possible to tag the initial flavor $`\alpha `$ and the final flavor $`\beta `$ of the neutrino. The flavor-mixing amplitude $`𝒜_{\alpha \beta }(T,𝐋)`$ is defined as the transition amplitude corresponding to the propagation over a time $`T`$ and a distance $`𝐋`$ of a neutrino initially of flavor $`\alpha `$, but detected with a flavor $`\beta `$.

The flavor-mixing amplitude can be expressed as a linear combination of amplitudes $`𝒜_j`$ corresponding to the propagation of different mass eigenstates:

$$𝒜_{\alpha \beta }(T,𝐋)=\underset{j}{}V_{\beta j}^{}𝒜_jV_{j\alpha },$$
(4)
where $`𝒜_j`$ is given by Eq. (2).

As the propagation time is not measured in experiments, the transition probability has to be averaged over $`T`$:

$$𝒫_{\alpha \beta }(𝐋)\underset{i,j}{}V_{i\alpha }V_{i\beta }^{}V_{j\alpha }^{}V_{j\beta }𝑑T𝒜_i𝒜_j^{}.$$
(5)
A discussion of the proportionality factor can be found in Ref. .

The next task consists in performing the momentum integration in the amplitude (2). This cannot be done without doing some assumptions on the overlap function (3), which should be as weak as possible. A fairly general oscillation formula can be obtained by approximating the in- and outgoing particles with gaussian wave packets . A gaussian wave packet is defined by

$$\psi _{\mathrm{SS}P_{in}}(𝐪,𝐐)=\left(\frac{2\pi }{\sigma _{p\mathrm{SS}P_{in}}^2}\right)^{3/4}\mathrm{exp}\left(\frac{(𝐪𝐐)^2}{4\sigma _{p\mathrm{SS}P_{in}}^2}\right),$$

where $`\sigma _{p\mathrm{SS}P_{in}}`$ is the width of the wave packet in momentum space. It is also useful to define the width $`\sigma _{x\mathrm{SS}P_{in}}`$ in configuration space by $`\sigma _{p\mathrm{SS}P_{in}}\sigma _{x\mathrm{SS}P_{in}}=1/2`$. The normalization of the states is defined as in Ref. , so that $`\frac{d𝐪}{(2\pi )^3}|\psi (𝐪)|^2=1`$.

If the wave packet is sharply peaked around its average momentum $`𝐐`$, i.e. $`\sigma _{\mathrm{SS}P_{in}}E_{\mathrm{SS}P_{in}}`$, the energy can be approximated by

$$E_{\mathrm{SS}P_{in}}(𝐪)E_{\mathrm{SS}P_{in}}(𝐐)+𝐯_{\mathrm{SS}P_{in}}(𝐪𝐐),$$

where $`E_{P_{in}}(𝐪)=\sqrt{𝐪^2+m_{\mathrm{SS}P_{in}}^2}`$ and $`𝐯_{\mathrm{SS}P_{in}}=𝐐/E_{\mathrm{SS}P_{in}}(𝐐)`$. This approximation means that dispersion is neglected in the external wave packets. The factors $`M_P(q,k)`$ and $`M_D(q^{},k^{})`$ multiplying the exponential vary slowly over the width of the wave packet and can be approximated by their value at the average momentum. They can be factorized outside the sum over the mass eigenstates since we have assumed that the neutrinos are either relativistic or nearly degenerate in mass.

With these approximations, the momentum integrations in the overlap function (3) can be done analytically. One then performs the integrations on $`𝐱`$, $`𝐱^{}`$ and $`t`$, $`t^{}`$, which are all gaussian. The final result can be written as

$$\psi (p^0,𝐩)=N\mathrm{exp}\left(f_{\mathrm{SS}P}(p^0,𝐩)f_{\mathrm{SS}D}(p^0,𝐩)\right),$$
(6)
with

$$f_{\mathrm{SS}P}(p^0,𝐩)=\frac{\left(𝐩𝐩_{\mathrm{SS}P}\right)^2}{4\sigma _{p\mathrm{SS}P}^2}+\frac{\left(p^0E_{\mathrm{SS}P}(𝐩𝐩_{\mathrm{SS}P})𝐯_{\mathrm{SS}P}\right)^2}{4\sigma _{e\mathrm{SS}P}^2},$$

where $`E_{\mathrm{SS}P}=E_{\mathrm{SS}P_{in}}E_{\mathrm{SS}P_{out}}`$, $`𝐩_{\mathrm{SS}P}=𝐐𝐊`$. The function $`f_{\mathrm{SS}D}(p^0,𝐩)`$ is defined in the same way, with the index $`P`$ replaced by $`D`$, except for the energy-momentum which is defined so as to be positive: $`E_{\mathrm{SS}D}=E_{\mathrm{SS}D_{out}}E_{\mathrm{SS}D_{in}}`$ and $`𝐩_{\mathrm{SS}D}=𝐊^{}𝐐^{}`$. The constant $`N`$ includes the normalization constants as well as the factors $`M_{P,D}`$ evaluated at the average momenta.

A new width $`\sigma _{p\mathrm{SS}P}`$ has been defined by $`\sigma _{p\mathrm{SS}P}\sigma _{x\mathrm{SS}P}=1/2`$, with

$$\frac{1}{\sigma _{x\mathrm{SS}P}^2}=\frac{1}{\sigma _{x\mathrm{SS}P_{in}}^2}+\frac{1}{\sigma _{x\mathrm{SS}P_{out}}^2}.$$

$`\sigma _{p\mathrm{SS}P}`$ can be interpreted as the momentum uncertainty at the source. The spatial width $`\sigma _{x\mathrm{SS}P}`$ is mainly determined by the external particle with the smallest space width. This is expected since the production region depends on the overlap in space-time of the external wave packets.

The symbol $`𝐯_{\mathrm{SS}P}`$ is defined by

$$𝐯_{\mathrm{SS}P}=\sigma _{x\mathrm{SS}P}^2\left(\frac{𝐯_{\mathrm{SS}P_{in}}}{\sigma _{x\mathrm{SS}P_{in}}^2}+\frac{𝐯_{\mathrm{SS}P_{out}}}{\sigma _{x\mathrm{SS}P_{out}}^2}\right).$$

It can be interpreted as the velocity of the production region, approximately equal to the velocity of the particle with the smallest spatial spread.

The symbol $`\mathrm{\Sigma }_{\mathrm{SS}P}`$, satisfying $`0\mathrm{\Sigma }_{\mathrm{SS}P}1`$, is defined by

$$\mathrm{\Sigma }_{\mathrm{SS}P}=\sigma _{x\mathrm{SS}P}^2\left(\frac{𝐯_{\mathrm{SS}P_{in}}^2}{\sigma _{x\mathrm{SS}P_{in}}^2}+\frac{𝐯_{\mathrm{SS}P_{out}}^2}{\sigma _{x\mathrm{SS}P_{out}}^2}\right).$$

Finally, the quantity

$$\sigma _{e\mathrm{SS}P}^2=\sigma _{p\mathrm{SS}P}^2\left(\mathrm{\Sigma }_{\mathrm{SS}P}𝐯_{\mathrm{SS}P}^2\right)\sigma _{p\mathrm{SS}P}^2$$

can be interpreted as the energy uncertainty at the source, or also as the inverse of the time of overlap of wave packets during the production process. Indeed, one can show that

$$\sigma _{e\mathrm{SS}P}^2=\underset{\alpha <\beta }{}\frac{\sigma _{x\mathrm{SS}P}^2}{4\sigma _{x\alpha }^2\sigma _{x\beta }^2}\left(𝐯_\alpha 𝐯_\beta \right)^2,$$
(7)
where the sum is over all wave packets involved in the production process. This sum is dominated by the term including the two smallest wave packets in configuration space (unless their velocities are nearly equal). If $`\sigma _{x1}`$ is the smallest width and $`\sigma _{x2}`$ the second smallest, one obtains

$$\sigma _{e\mathrm{SS}P}\frac{|𝐯_1𝐯_2|}{\sigma _{x2}}\frac{1}{T_{\mathrm{SS}P}^{overlap}},$$
(8)
where $`T_{\mathrm{SS}P}^{overlap}`$ is defined as the duration of the production process. Thus, $`\sigma _{e\mathrm{SS}P}`$ can be interpreted as the energy uncertainty at the source, since it is proportional to the inverse of the time of overlap of the external wave packets at the source. The quantities $`\sigma _{x\mathrm{SS}D}`$, $`\sigma _{p\mathrm{SS}D}`$, $`𝐯_{\mathrm{SS}D}`$, $`\mathrm{\Sigma }_{\mathrm{SS}D}`$, $`\sigma _{e\mathrm{SS}D}`$, $`T_{\mathrm{SS}D}^{overlap}`$ have similar definitions and properties.

Note that stationary boundary conditions are recovered by setting $`𝐯_{\mathrm{SS}P,D}=0`$ and $`\sigma _{e\mathrm{SS}P,D}=0`$, with $`\sigma _{p\mathrm{SS}P,D}`$ different from zero. Moreover, it is reasonable to impose the constraints $`|𝐯_{\mathrm{SS}P,D}|T_{\mathrm{SS}P,D}^{overlap}S_{\mathrm{SS}P,D}`$, where $`S_{\mathrm{SS}P}`$ (resp. $`S_{\mathrm{SS}D}`$) is the size of the macroscopic region of production (resp. detection). Thus we shall assume the following constraints:

$$\frac{|𝐯_{\mathrm{SS}P,D}|}{\sigma _{e\mathrm{SS}P,D}}S_{\mathrm{SS}P,D}.$$
(9)
These bounds are very conservative, since we shall see that stationary models such as those found in Refs.  are recovered by setting $`|𝐯_{\mathrm{SS}P,D}|/\sigma _{e\mathrm{SS}P,D}=0`$. In the example of the Grimus-Stockinger model , an initial stationary neutron ($`𝐯_{\mathrm{SS}P_{in},n}=0`$) decays into a stationary proton ($`𝐯_{\mathrm{SS}P_{out},pr}=0`$), a ‘plane-wave’ electron ($`\sigma _{x\mathrm{SS}P_{out},el}=\mathrm{}`$) and the intermediate neutrino. At detection, the neutrino collides with a stationary electron ($`𝐯_{\mathrm{SS}D_{in}}=0`$) and the outgoing neutrino and electron are represented as plane waves ($`\sigma _{x\mathrm{SS}D_{out},\nu }=\sigma _{x\mathrm{SS}D_{out},el}=\mathrm{}`$).

As the propagation distance is macroscopic, only processes satisfying global conservation of energy-momentum have a non-negligible probability of occurring. Since our aim is not to prove this well-known fact, we impose that

$$𝐩_{\mathrm{SS}P}=𝐩_{\mathrm{SS}D}𝐩_0\text{and}E_{\mathrm{SS}P}=E_{\mathrm{SS}D}E_0.$$
(10)
This approximation allows to do expansions around $`𝐩_0`$ and $`E_0`$. An associated velocity can be defined by $`𝐯_0=𝐩_0/E_0`$.

## 3 The oscillation formula, with or without stationary boundary conditions

In this section, we shall explain a quick method to evaluate the flavor-changing probability (5) with the amplitude given by Eq. (2) and the overlap function by Eq. (6). Since the experimental conditions are such that the propagating particle is on-shell, the main contribution to the transition amplitude (2) comes from the pole of the propagator. However, one has to be careful with the choice of the contour as the analytic continuation of most overlap functions diverges at infinity in the complex plane. The integration on the 3-momentum can be done with the help of the Grimus-Stockinger theorem . Let $`\psi (𝐩)`$ be a 3 times continuously differentiable function on $`𝐑^3`$ such that $`\psi `$ itself and all its first and second derivatives decrease at least like $`1/𝐩^2`$ for $`|𝐩|\mathrm{}`$. Then, for a real number $`A>0`$,

$$d^3p\frac{\psi (𝐩)e^{i𝐩𝐋}}{A𝐩^2+iϵ}\stackrel{L\mathrm{}}{}\frac{2\pi ^2}{L}\psi (\sqrt{A}𝐥)e^{i\sqrt{A}L},$$

where $`L=|𝐋|`$ and $`𝐥=𝐋/L`$. For $`A<0`$, the integral decreases like $`L^2`$.

The remaining energy integral in the amplitude (2) can be done by a saddle-point approximation . However, it is quicker to perform first the time average in the probability (5), which yields a delta function, and makes one of the energy integrations trivial. At this stage, one has

$$𝑑T𝒜_i𝒜_j^{}\frac{1}{L^2}𝑑E\psi (E,q_i𝐥)\psi ^{}(E,q_j𝐥)e^{i(q_iq_j)L},$$
(11)
where $`\psi (E,𝐩)`$ is the overlap function defined by Eq. (6) and $`q_j=\sqrt{E^2m_j^2}`$.

Eq. (11) shows that the transition probability can be interpreted as an incoherent sum (i.e. occurring in the probability) over energy eigenstates: interference occurs only between the components of $`\psi (E,𝐩)`$ having the same energy . In this way, the correspondence between models with and without stationary boundary conditions is obvious: the time-integrated nonstationary probability is equivalent to the energy-integrated stationary probability. For example, the oscillation formula obtained by Grimus and Stockinger with stationary boundary conditions  has the form of the integrand in the right-hand side of Eq. (11). This equivalence confirms that the stationary case can be obtained from the more general nonstationary case in the limit of a vanishing energy width. This limit will have to be checked explicitly on the final oscillation formula, as it is not obvious that the approximations involved in the computations respect this feature.

If the coordinate system is chosen so that $`𝐋`$ is oriented along a coordinate axis, it is easy to rewrite the integral (11) as

$$𝑑T𝒜_i𝒜_j^{}\frac{g^2(𝐥)}{L^2}𝑑Ee^{i(q_iq_j)Lf_i(E)f_j(E)},$$
(12)
with the definitions $`f_j(E)=f_{j\mathrm{SS}P}(E)+f_{j\mathrm{SS}D}(E)`$ and

$$f_{j\mathrm{SS}P}(E)=\frac{\left(\sqrt{E^2m_j^2}p_0\right)^2}{4\sigma _{p\mathrm{SS}P}^2}+\frac{\left(EE_0\left(\sqrt{E^2m_j^2}p_0\right)v_{\mathrm{SS}P}\right)^2}{4\sigma _{e\mathrm{SS}P}^2},$$
(13)
where $`v_{\mathrm{SS}P}`$ and $`p_0`$ are the components of $`𝐯_{\mathrm{SS}P}`$ and $`𝐩_0`$ along $`𝐋`$, while $`E_0`$ has been redefined so as to absorb the transversal part of $`𝐩_0𝐯_{\mathrm{SS}P}`$. The definition of $`f_{j\mathrm{SS}D}(E)`$ is similar. The function $`g(𝐥)`$ expresses the geometrical constraint between the direction of observation $`𝐥=𝐋/L`$ and the momentum $`𝐩_0`$:

$$g(𝐥)=\mathrm{exp}\left(\frac{(𝐩_0\times 𝐥)^2}{4\sigma _p^2}\right),$$
(14)
and restricts the neutrino propagation to a cone of axis $`𝐩_0`$ and angle $`arcsin(\sigma _p/p_0)`$. The momentum width $`\sigma _p`$, defined by

$$\frac{1}{\sigma _p^2}=\frac{1}{\sigma _{p\mathrm{SS}P}^2}+\frac{1}{\sigma _{p\mathrm{SS}D}^2},$$
(15)
is approximately equal to the smallest width among the production and detection momentum widths. The associated width $`\sigma _x`$ in configuration space is defined by $`\sigma _p\sigma _x=1/2`$.

The remaining energy integral in Eq. (12) can be performed as a gaussian integral by expanding the integrand to second order around the maximum of its modulus (this is called Laplace’s method ). The result takes a much simpler form if an expansion in small mass differences is done around ($`E_0`$, $`p_0`$, $`m_0`$, $`v_0`$), where $`m_0^2=E_0^2p_0^2`$ and $`v_0=p_0/E_0`$. The expansion parameter is noted $`ϵ`$ and refers collectively to all $`\delta m_j^2/2E_0^2`$, where $`\delta m_j^2=m_j^2m_0^2`$. We shall calculate the transition probability to $`𝒪(ϵ^2)`$ in the real part of the argument of the exponential (since the order $`𝒪(ϵ)`$ vanishes) and to $`𝒪(ϵ)`$ in the phase. The gaussian integration will be consistent with the $`ϵ`$-expansion if the extremum is computed to $`𝒪(ϵ)`$, the real part of the argument of the exponential to $`𝒪(ϵ^2)`$, the phase to $`𝒪(ϵ)`$, the first derivatives to $`𝒪(ϵ)`$ and the second derivatives to $`𝒪(ϵ^0)`$. It is important to take into account all these terms, if the result is to coincide with the results obtained by the methods explained in the fourth section.

The modulus of the integrand in Eq. (12) is maximal for

$$E_{ij}=E_0+\rho \frac{\delta m_i^2+\delta m_j^2}{4E_0}+𝒪(ϵ^2),$$

where the dimensionless number $`\rho `$ is defined by

$$\rho =\sigma _{peff}^2\left(\frac{1}{\sigma _p^2}\frac{v_{\mathrm{SS}P}(v_0v_{\mathrm{SS}P})}{\sigma _{e\mathrm{SS}P}^2}\frac{v_{\mathrm{SS}D}(v_0v_{\mathrm{SS}D})}{\sigma _{e\mathrm{SS}D}^2}\right).$$
(16)
The effective width in momentum space $`\sigma _{peff}`$ is defined by

$$\frac{1}{\sigma _{peff}^2}=\frac{1}{\sigma _p^2}+\frac{(v_0v_{\mathrm{SS}P})^2}{\sigma _{e\mathrm{SS}P}^2}+\frac{(v_0v_{\mathrm{SS}D})^2}{\sigma _{e\mathrm{SS}D}^2},$$
(17)
with $`\sigma _p`$ defined by Eq. (15). The effective width can be interpreted as the energy-momentum width of the oscillating particle, since it is the width of the overlap function. Indeed, the value to $`𝒪(ϵ^0)`$ of the second derivative of $`f_i(E)+f_j(E)`$ at the extremum reads

$$\frac{1}{2}\frac{d^2(f_i+f_j)}{dE^2}(E_{ij})=\frac{1}{2v_0^2\sigma _{peff}^2}.$$
(18)
The effective width is dominated by the smallest among the energy uncertainties (recall that $`\sigma _{e\mathrm{SS}P,D}\sigma _{p\mathrm{SS}P,D}`$). The effective width $`\sigma _{xeff}`$ in configuration space, defined by $`\sigma _{peff}\sigma _{xeff}=1/2`$, is then approximately equal either to the production or to the detection time uncertainty, depending on which one is the largest. In the stationary limit ($`𝐯_{\mathrm{SS}P,D}=0`$ and $`\sigma _{e\mathrm{SS}P,D}=0`$), $`\sigma _{peff}`$ goes to zero.

The parameter $`\rho `$ has been defined so as to be in correspondence with the notation of Ref. . The symbol $`\omega `$ appearing in that article is related to our notation by $`\omega =\sigma _p^2/\sigma _{peff}^2`$. Note that the authors of Ref.  do not compute $`\rho `$ explicitly and also take the relativistic limit $`v_0=1`$. The explicit value of $`\rho `$ is very interesting to know, since $`\rho =0`$ in the case of stationary boundary conditions, in which case all mass eigenstates have the same energy $`E_0`$.

The value to $`𝒪(ϵ^2)`$ of $`f_i(E)+f_j(E)`$ at the extremum reads

$$f_i(E_{ij})+f_j(E_{ij})=\frac{(\delta m_i^2)^2+(\delta m_j^2)^2}{16\sigma _m^2E_0^2}+2\pi ^2\left(\frac{\rho \sigma _{xeff}}{L_{ij}^{osc}}\right)^2,$$
(19)
where the mass width $`\sigma _m`$ is defined by

$$\frac{1}{\sigma _m^2}=\sigma _{peff}^2\left(\frac{1}{\sigma _p^2}\left(\frac{1}{\sigma _{e\mathrm{SS}P}^2}+\frac{1}{\sigma _{e\mathrm{SS}D}^2}\right)+\frac{(v_{\mathrm{SS}P}v_{\mathrm{SS}D})^2}{\sigma _{e\mathrm{SS}P}^2\sigma _{e\mathrm{SS}D}^2}\right).$$
(20)
The expansion of the phase present in Eq. (12) around the extremum reads

$$(q_iq_j)L\frac{\delta m_{ij}^2L}{2p_0}+\frac{\delta m_{ij}^2L}{2p_0^2v_0}(EE_{ij}).$$
(21)
The second derivative of the phase is of $`𝒪(ϵ)`$ and can be neglected with respect to the second derivative of $`f_i+f_j`$.

The gaussian integration in Eq. (12) may now be performed with the help of Eqs. (18), (19) and (21). As a result, the flavor-mixing transition probability (5) corresponding to the propagation over a distance $`𝐋=L𝐥`$ of a neutrino, initially of flavor $`\alpha `$, but detected with a flavor $`\beta `$, reads

$`𝒫_{\alpha \beta }(L𝐥){\displaystyle \frac{1}{L^2}}\mathrm{exp}\left({\displaystyle \frac{(𝐩_0\times 𝐥)^2}{2\sigma _p^2}}\right){\displaystyle \underset{i,j}{}}V_{i\alpha }V_{i\beta }^{}V_{j\alpha }^{}V_{j\beta }`$
$`\times \mathrm{exp}\left(2\pi i{\displaystyle \frac{L}{L_{ij}^{osc}}}\left({\displaystyle \frac{L}{L_{ij}^{coh}}}\right)^22\pi ^2\left({\displaystyle \frac{\rho \sigma _{xeff}}{L_{ij}^{osc}}}\right)^2{\displaystyle \frac{(\delta m_i^2)^2+(\delta m_j^2)^2}{16\sigma _m^2E_0^2}}\right).`$ (22)
The oscillation length $`L_{ij}^{osc}`$ for the masses $`m_i`$ and $`m_j`$ is given by

$$L_{ij}^{osc}=\frac{4\pi p_0}{\delta m_{ij}^2},$$
(23)
where $`\delta m_{ij}^2=m_i^2m_j^2`$ is taken to be positive. The coherence length $`L_{ij}^{coh}`$ is defined by

$$L_{ij}^{coh}=\frac{1}{\sqrt{2}\pi }\frac{p_0}{\sigma _{peff}}L_{ij}^{osc}.$$
(24)
Let us now proceed to the analysis of the oscillation formula (22). The first term in the second exponential of Eq. (22) is the standard oscillation phase proportional to the propagation distance. The second term, or coherence-length term, leads to the vanishing of oscillations beyond the coherence length $`L_{ij}^{coh}`$. This phenomenon is partly due to the progressive separation of mass-eigenstates wave packets propagating in space (see discussion after Eq. (68)). The third term is a localization term, i.e. a constraint imposing that oscillations vanish unless the oscillation length is larger than the space-time uncertainty:

$$L_{ij}^{osc}\rho \sigma _{xeff}.$$
(25)
This condition can be rewritten as $`\delta m_{ij}^2/p_0\sigma _{peff}/\rho `$, stating that oscillations vanish if the energy-momentum measurements allow to distinguish between the different mass eigenstates. The coherence length and localization conditions were already predicted in intermediate wave packet models of neutrino oscillations .

The fourth term in the second exponential of Eq. (22) could be a matter of concern since it does not vanish in the limit $`m_i=m_j`$. First of all, note that this kind of term is not specific to our computation. For example, it would appear in the oscillation formula (26) of Ref.  if the terms $`S_a(E_a)`$ present in Eq. (22) of that article had been expanded beyond zeroth order in $`m_a^2/E_a^2`$. It is of interest to rewrite the term under discussion as

$$\frac{(\delta m_i^2)^2+(\delta m_j^2)^2}{16\sigma _m^2E_0^2}=\frac{(\delta m_{ij}^2)^2}{32\sigma _m^2E_0^2}+\frac{(\delta m_i^2+\delta m_j^2)^2}{32\sigma _m^2E_0^2}.$$
(26)
The first term in Eq. (26) is recognized as a localization constraint and can be rewritten as

$$L_{ij}^{osc}\sigma _x,$$
(27)
as $`\sigma _mv_0\sigma _p`$ whether the stationary limit is taken or not. The second term in Eq. (26) imposes that

$$\frac{|m_i^2+m_j^22m_0^2|}{E_0}\sigma _m.$$
(28)
As $`m_0`$ is related to the average in- and outgoing momentum by $`m_0^2=E_0^2p_0^2`$, condition (28) means that the mass eigenstates have to be on-shell with respect to $`(E_0,p_0)`$ within the uncertainty $`\sigma _m`$. For example, this constraint is impossible to satisfy in the mixing of relativistic and nonrelativistic neutrinos if the available energy-momentum is such that only the lightest neutrino can be produced. However it has no effect on the oscillations in the two cases considered in the present article, namely relativistic neutrinos or nearly mass degenerate neutrinos. Condition (28) should simply be considered as expressing the conservation of energy-momentum. Such kinematical constraints are usually not included in the oscillation formula, though they rightly belong to it. A complete computation of the transition probability should not only include this energy-momentum constraint, but also similar terms arising from the interaction amplitudes $`M_{P,D}`$, from the prefactor resulting from the gaussian integration and from the spin structure of the propagator. Actually the neutrino masses should be expected to appear not only through mass differences but also through their absolute values.

The oscillation formula (22) is similar to Eq. (26) of Ref. . However, three interesting new elements appear in Eq. (22): the parameter $`\rho `$ has been computed explicitly, the oscillation formula (22) is valid for relativistic or nonrelativistic (but nearly mass degenerate) neutrinos, and the dispersion of the ‘wave packet’ corresponding to the neutrino has been taken into account. This last fact is not at all obvious from the above computation, in which the amplitude $`𝒜_j`$ has been inserted into the probability before doing the integration over the energy. Had this energy integration been performed first, the explicit dependence of $`𝒜_j`$ on $`(T,𝐋)`$ would have been obtained, i.e. a gaussian spreading with $`T`$ in the direction $`𝐋`$. In section 4.3 this spreading will be interpreted as the longitudinal dispersion of the ‘wave packet’ associated with the propagating neutrino. We shall also see that computations neglecting this dispersion lead to a result differing from Eq. (22) in the nonrelativistic limit (see Eq. (67)). In contrast to the method of section 4.3, the question of dispersion does not even arise in the derivation of Eqs. (11)-(22) in which it is automatically taken into account.

Another feature of the oscillation formula (22) is the increase of the coherence length when a long coherent measurement in time is performed at the detector, even if the neutrino ‘wave packets’ have separated spatially . In that case, the energy uncertainty at detection goes to zero, $`\sigma _{e\mathrm{SS}D}0`$, so that the effective width also goes to zero, $`\sigma _{peff}0`$, and the coherence length becomes infinite, $`L_{ij}^{coh}\mathrm{}`$. At first sight, oscillations seem to vanish in that limit, because the localization term depending on $`\rho \sigma _{xeff}`$ seems to diverge . If it were true, it would be impossible to increase without limit the coherence length by performing long coherent measurements. Note that this would be in contradiction with stationary boundary condition models, which have a zero energy uncertainty but an infinite coherence length. This apparent contradiction can be cleared up by examining carefully the term $`\rho \sigma _{xeff}`$ in the limit $`\sigma _{e\mathrm{SS}D}0`$. With the help of Eqs. (9), (16) and (17), one sees that

$$|\rho |\sigma _{xeff}\stackrel{\sigma _{e\mathrm{SS}D}0}{}\frac{|v_{\mathrm{SS}D}|}{\sigma _{e\mathrm{SS}D}}S_{\mathrm{SS}D},$$

where $`S_{\mathrm{SS}D}`$ was defined as the size of the macroscopic detection region. Thus the localization term does not give a stronger constraint than $`S_{\mathrm{SS}D}L_{ij}^{osc}`$. This constraint is always satisfied, as it is equivalent to the constraint obtained by averaging the transition probability over the production region. Therefore, the coherence length can be increased without bound by more accurate energy measurements, contrary to what was claimed in Refs. . Note that this is not true if the accuracy of the 3-momentum measurements is increased, as the localization makes the oscillations vanish when the corresponding spatial uncertainty becomes larger than the oscillation length. Of course the opposite conclusions would be reached if experiments measured time, not distance.

As a corollary, the stationary limit of the oscillation formula (22) is well-defined. Stationary boundary conditions, $`v_{\mathrm{SS}P,D}=0`$ and $`\sigma _{e\mathrm{SS}P,D}=0`$, lead to an infinite effective width $`\sigma _{xeff}`$ and thus to an infinite coherence length. However, the product $`\rho \sigma _{xeff}`$ remains finite, as we have just explained. For example, the Grimus-Stockinger model  is obtained in the limit $`\rho =0`$, $`\sigma _{xeff}\mathrm{}`$ with $`\rho \sigma _{xeff}=0`$. The latter condition means that this model can be recovered from the external wave packet model if the stationary limit has the property $`v_{\mathrm{SS}P,D}/\sigma _{e\mathrm{SS}P,D}=0`$ (see discussion after Eq. (9)). Though condition (25) becomes ineffective in that limit, the localization condition (27) is still present. Thus there is no contradiction between models with stationary boundary conditions and those with external wave packets. The former type of model can be obtained from the latter in a smooth limit.

In conclusion, if the observability conditions $`|𝐩_0\times 𝐥|\sigma _p`$, $`LL_{ij}^{coh}`$ and $`L_{ij}^{osc}\sigma _x`$ are satisfied, the oscillation formula (22) reduces to the standard formula (1), with the additional property of $`1/L^2`$ geometrical decrease.

## 4 Other oscillation formulas

### 4.1 The Jacob-Sachs theorem

With result (22), we have been able to reconcile the oscillation formulas proposed on the one hand by Giunti, Kim and Lee  and on the other hand by Grimus and Stockinger . More generally, the stationary boundary conditions models have been shown to be special cases of the external wave packet models. We have seen that Eq. (22) reduces to the standard oscillation formula (1) if some observability conditions are verified; otherwise oscillations vanish. However, other authors also working in a quantum field framework have argued that formula (22) is not valid for all experimental conditions.

First, Ioannisian and Pilaftsis  derive a formula for neutrino oscillations which exhibits a plane wave behaviour, if the condition $`L/\sigma _xp_0/\sigma _p`$ is satisfied (with $`\sigma _xL_{ij}^{osc}`$, as usual). By ‘plane wave oscillations’, these authors mean that the oscillation amplitude $`𝒜_j`$ does not decrease as the inverse of the distance, that its phase depends on the direction $`𝐋`$ as $`𝐩_𝐣𝐋`$, and that the amplitude is not negligible in directions other than along the neutrino momentum. If it were true, the oscillation length measured by a specific detector would depend on the direction of the total momentum of the initial particles. To give a typical example, the above ‘plane wave condition’ is satisfied for atmospheric neutrinos if $`\sigma _x10^2`$cm (with $`L1000`$km and $`p1`$GeV). Unfortunately, this condition appears nowhere in the treatment of the previous section. Actually, it is not clear whether the Grimus-Stockinger theorem is valid under this condition since its derivation uses the large $`L`$ limit.

Secondly, Shtanov  argues that the relative weights of the oscillating exponentials can be strongly sensitive to the neutrino masses, if the source and detector are strongly localized. Since this author works in configuration space, the derivation of the previous section cannot be useful to assess his claim.

Thus, the examination of the claims made by Ioannisian, Pilaftsis and Shtanov requires another method of evaluation of the amplitude (2). Instead of using the Grimus-Stockinger theorem to integrate on the 3-momentum, we shall first perform the integration on the energy with the Jacob-Sachs theorem . This theorem is based on the assumption that the energy spectrum of all incident particles is limited to a finite range. Thus the overlap function $`\psi (E,𝐩)`$ is distinct from zero only for $`p^2=E^2𝐩^2`$ within certain bounds (with $`E>0`$). On this interval, $`\psi (E,𝐩)`$ is taken to be infinitely differentiable. In that case, the Jacob-Sachs theorem says that the asymptotic value of the energy integral in Eq. (2), when $`T\mathrm{}`$, is given by its residue at the pole below the real axis. Thus the evaluation of the partial amplitude (2) with this theorem yields

$$𝒜_j\frac{\pi }{(2\pi )^4}\frac{d^3p}{\sqrt{𝐩^2+m_j^2}}\psi (\sqrt{𝐩^2+m_j^2},𝐩)e^{i\varphi _j(𝐩)},$$
(29)
where $`\varphi _j(𝐩)=\sqrt{𝐩^2+m_j^2}T𝐩𝐋`$. For gaussian external wave packets, the overlap function is given by Eq. (6). In principle, this function should be cut off outside the energy range determined by experimental conditions so as to satisfy the conditions of the Jacob-Sachs theorem. However these corrections are very small and will be neglected in the computations.

It is interesting to note that the amplitude (29) is mathematically equivalent to the amplitude obtained in the intermediate wave packet model , in which the neutrino mass eigenstates are directly represented by wave packets. The overlap function $`\psi (\sqrt{𝐩^2+m_j^2},𝐩)`$ corresponds to the wave function of the jth mass eigenstate. Thus, it makes sense, in an external wave packet model, to talk about mass eigenstate ‘wave packets’ associated with the propagating neutrino. Note however that the overlap function takes into account not only the properties of the source, but also of the detector, which is a bit unusual for a wave packet interpretation.

### 4.2 No-dispersion regime

#### 4.2.1 Choice of integration method

The integration over the 3-momentum in Eq. (29) cannot be done analytically. Resorting to the explicit form (6) of the overlap function valid for gaussian external wave packets, we see that the integral (29) can be approximated by means of an asymptotic expansion for which two kinds of large parameters can be used. On the one hand, $`\sigma _{p\mathrm{SS}P,D}^2`$ and $`\sigma _{e\mathrm{SS}P,D}^2`$ are large parameters appearing in the overlap function (6). They suggest a second order expansion of the integrand around the maximum $`𝐩_j`$ of the overlap function, followed by a gaussian integration: this is called Laplace’s method . On the other hand, $`T`$ and $`𝐋`$ are large parameters appearing in the phase. They suggest a second order expansion of the integrand around the stationary point $`𝐩_{cl,j}`$ of the phase, followed by a gaussian integration: this is called the method of stationary phase . The competition between these two asymptotic behaviors implies a detailed study of the oscillation of the phase around the average momentum $`𝐩_j`$. The expansion of the phase in Eq. (29) should be compared with the expansion of the overlap function. Although both methods are expected to lead roughly to the same answer in the case of the propagation of a single particle, it should be checked whether the delicate compensation mechanism resulting in the oscillation phase is independent of the method chosen.

The study of the overlap function amounts to the study of the argument of the exponential, i.e. of the function $`(f_{\mathrm{SS}P}+f_{\mathrm{SS}D})(𝐩)`$ (see Eq. (6) with $`p^0=\sqrt{𝐩^2+m_j^2}`$). Using as in the previous section an expansion in small mass differences $`ϵ=\delta m_j^2/E_0^2`$, the value $`𝐩_j`$ minimizing $`f_{\mathrm{SS}P}+f_{\mathrm{SS}D}`$ is given to $`𝒪(ϵ)`$ by

$$𝐩_j=𝐩_0+\left(\alpha 𝐮_{\mathrm{SS}P}+\beta 𝐮_{\mathrm{SS}D}\right)\frac{\delta m_j^2}{2E_0},$$
(30)
where

$$𝐮_{\mathrm{SS}P,D}=\frac{𝐯_0𝐯_{\mathrm{SS}P,D}}{2\sigma _{e\mathrm{SS}P,D}}.$$

The associated energy $`E_j=\sqrt{𝐩_j^2+m_j^2}`$ can be expanded to $`𝒪(ϵ)`$ as

$$E_j=E_0+\stackrel{~}{\rho }\frac{\delta m_j^2}{2E_0},$$
(31)
where $`\stackrel{~}{\rho }=1+\alpha 𝐯_0𝐮_{\mathrm{SS}P}+\beta 𝐯_0𝐮_{\mathrm{SS}D}`$. The values of the dimensionless coefficients $`\alpha `$ and $`\beta `$ can be computed but we shall not need their explicit expressions. It is sufficient to know that $`\stackrel{~}{\rho }0`$ in the stationary limit. A velocity $`𝐯_j=𝐩_j/E_j`$ is also defined for future use.

We are now going to approximate the overlap function as a gaussian and compute its three characteristic widths. At the extremum $`𝐩_j`$, the Hessian matrix of $`f_{\mathrm{SS}P}+f_{\mathrm{SS}D}`$ reads to $`𝒪(ϵ^0)`$

$`\mathrm{\Sigma }^{ab}`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2(f_{\mathrm{SS}P}+f_{\mathrm{SS}D})}{p^ap^b}}(𝐩_j)`$
$`=`$ $`{\displaystyle \frac{\delta ^{ab}}{4\sigma _p^2}}+u_{\mathrm{SS}P}^au_{\mathrm{SS}P}^b+u_{\mathrm{SS}D}^au_{\mathrm{SS}D}^b,`$
where $`\sigma _p`$ is defined by Eq. (15). The matrix $`\mathrm{\Sigma }^{ab}`$ determines the range of $`𝐩`$ values for which the overlap function $`\psi (\sqrt{𝐩^2+m_j^2},𝐩)`$ is not negligible. As $`\mathrm{\Sigma }^{ab}`$ is symmetric, it can be diagonalized by an orthogonal coordinate transformation. The eigenvalues of $`\mathrm{\Sigma }^{ab}`$ are

$`\sigma _x^2`$ $`=`$ $`{\displaystyle \frac{1}{4\sigma _p^2}},`$
$`\sigma _{x\pm }^2`$ $`=`$ $`{\displaystyle \frac{1}{4\sigma _p^2}}+{\displaystyle \frac{1}{2}}\left(𝐮_{\mathrm{SS}P}^2+𝐮_{\mathrm{SS}D}^2\right)\pm {\displaystyle \frac{1}{2}}\sqrt{\left(𝐮_{\mathrm{SS}P}^2+𝐮_{\mathrm{SS}D}^2\right)^24\left(𝐮_{\mathrm{SS}P}\times 𝐮_{\mathrm{SS}D}\right)^2}.`$
The eigenvector associated with $`\sigma _x^2`$ is in the direction of $`𝐮_{\mathrm{SS}P}\times 𝐮_{\mathrm{SS}D}`$, while the eigenvectors associated with $`\sigma _{x\pm }^2`$ belong to the plane defined by $`𝐮_{\mathrm{SS}P}`$ and $`𝐮_{\mathrm{SS}D}`$. In the limit $`|𝐮_{\mathrm{SS}P}||𝐮_{\mathrm{SS}D}|`$ (resp. $`|𝐮_{\mathrm{SS}P}||𝐮_{\mathrm{SS}D}|`$), the eigenvalues $`\sigma _x^2`$ and $`\sigma _x^2`$ become degenerate and the eigenvector associated with $`\sigma _{x+}^2`$ becomes aligned with $`𝐮_{\mathrm{SS}P}`$ (resp. $`𝐮_{\mathrm{SS}D}`$). This is also the case in the limit of parallel $`𝐮_{\mathrm{SS}P}`$ and $`𝐮_{\mathrm{SS}D}`$. These limits are relevant to the case of stationary boundary conditions which are examined below.

Let us choose coordinate axes $`(𝐞_x,𝐞_y,𝐞_z)`$ coinciding with the normalized eigenvectors associated with $`(\sigma _x^2,\sigma _x^2,\sigma _{x+}^2)`$ respectively. The quantities $`(\sigma _p^2,\sigma _p^2,\sigma _{p+}^2)`$ (with $`\sigma _{p\pm }\sigma _{x\pm }=1/2`$) can be interpreted as the momentum widths of the overlap function, since they give constraints on the range of $`𝐩`$ values for which the overlap function is non-negligible:

$`|p^xp_j^x|`$ $``$ $`\sigma _p,`$
$`|p^yp_j^y|`$ $``$ $`\sigma _p,`$
$`|p^zp_j^z|`$ $``$ $`\sigma _{p+}.`$ (32)
The case of the stationary limit is of special interest. Recall that stationary boundary conditions are obtained in the external wave packet model by taking $`𝐯_{\mathrm{SS}P,D}0`$ and $`\sigma _{e\mathrm{SS}P,D}0`$ with $`|𝐯_{\mathrm{SS}P,D}|\sigma _{e\mathrm{SS}P,D}S_{\mathrm{SS}P,D}`$ (see Eq. (9)). In this limit, the axis $`𝐞_z`$ becomes aligned with $`𝐯_0`$,

$$v_0^{x,y}v_{\mathrm{SS}P,D}^{x,y}0,$$
(33)
and two eigenvalues become degenerate while the third diverges:

$`\sigma _x^2`$ $``$ $`\sigma _x^2,`$
$`\sigma _{x+}^2`$ $``$ $`{\displaystyle \frac{1}{4\sigma _p^2}}+𝐮_{\mathrm{SS}P}^2+𝐮_{\mathrm{SS}D}^2,`$ (34)
In other words, the transversal widths (i.e. in the directions orthogonal to $`𝐩_0`$) are given in the stationary limit by $`\sigma _p`$, while the longitudinal width (i.e. in the direction of $`𝐩_0`$) is given by $`\sigma _{p+}\sigma _p`$. Note that the asymptotic value of $`\sigma _{x+}^2`$ is very similar to the definition (17) of $`\sigma _{xeff}^2`$.

The expansion of the overlap function has to be compared with the expansion of the phase $`\varphi _j(𝐩)`$ around $`𝐩_j`$, which reads

$$\varphi _j(𝐩)\varphi _j(𝐩_j)+(𝐯_jT𝐋)(𝐩𝐩_j)+\frac{T}{2E_0}(p^ap_j^a)R^{ab}(p^bp_j^b),$$
(35)
where $`R^{ab}=\delta ^{ab}v_0^av_0^b`$. The second derivatives have been evaluated to $`𝒪(ϵ^0)`$.

Laplace’s method will be preferred to the method of stationary phase if the phase $`\varphi _j(𝐩)`$ varies slowly over the ‘bump’ of the overlap function. As $`\sigma _p\sigma _p\sigma _{p+}`$, the variation of the phase will become important in the first place along the direction $`p^x`$, then in the direction $`p^y`$ and finally in the direction $`p^z`$. The criterion for the use of Laplace’s method in all three directions $`p^{x,y,z}`$ will thus be determined by considering the largest momentum width $`\sigma _p`$. The insertion of the constraints (32) into the phase (35) yields first order conditions for a slowly varying phase,

$`|v_j^xTL^x|\sigma _p`$ $``$ $`1,`$
$`|v_j^yTL^y|\sigma _p`$ $``$ $`1,`$
$`|v_j^zTL^z|\sigma _{p+}`$ $``$ $`1,`$ (36)
as well as a second order condition,

$$\frac{T}{E_0}\sigma _p^21,$$
(37)
where we have used the property $`q^aR^{ab}q^b𝐪^2`$ and the fact that $`\sigma _p`$ is the largest width. For a given $`T`$ satisfying Eq. (37), it is always possible to find a range of $`𝐋`$ values so that conditions (36) are satisfied. For other $`𝐋`$ values, the amplitude is negligible as will be checked on the result (see Eq. (42)). Thus the criterion allowing to choose the integration method is determined by Eq. (37): the integration on $`𝐩`$ will be done by Laplace’s method if $`TE_0/\sigma _p^2`$ or, equivalently (with $`L=|𝐋|`$ and $`p_0=|𝐩_0|`$), if

$$L\frac{p_0}{\sigma _p^2},$$
(38)
since conditions (36) impose the relation $`𝐋𝐯_0T`$ as long as $`\sigma _{x+}L`$.

For $`\sigma _{x+}L`$ (stationary limit), we now show that condition (38) is directly obtained without going through condition (37). The overlap function imposes in that limit that $`|𝐩|=\sqrt{E_0^2m_j^2}`$, so that we are left with an angular integration with the angular part of the integrand given by

$$\mathrm{exp}\left(\frac{𝐩𝐩_0}{2\sigma _p^2}+i𝐩𝐋\right).$$
(39)
Condition (38) shows that the angular variation of the phase in (39) is slow with respect to the angular variation of the overlap function, in which case Laplace’s method will give good results. Therefore condition (38) constitutes a good criterion for the use of Laplace’s method whether the stationary limit is taken or not.

Condition (38) is usually not verified in oscillation experiments, because $`L/\sigma _xp_0/\sigma _p`$ in most cases. This condition is the same than the one under which Ioannisian and Pilaftsis  obtain ‘plane wave’ oscillations.

#### 4.2.2 Amplitude

If condition (38) is satisfied, the evaluation of the integral (29) can be done by Laplace’s method and yields

$$𝒜_j\sigma _p\sigma _p\sigma _{p+}\mathrm{exp}\left(iE_jT+i𝐩_j𝐋\left(\frac{\delta m_j^2}{4\stackrel{~}{\sigma }_mE_0}\right)^2F_j(T)\right),$$
(40)
where the third term in the exponential comes from $`(f_{\mathrm{SS}P}+f_{\mathrm{SS}D})(𝐩_j)`$; the parameter $`\stackrel{~}{\sigma }_m`$ can be computed and has the dimension of a momentum width. The function $`F_j(T)`$ is defined by

$$F_j(T)=\frac{1}{4}(𝐯_jT𝐋)^t\left(\mathrm{\Sigma }+i\frac{T}{2E_0}R\right)^1(𝐯_jT𝐋),$$
(41)
where $`\mathrm{\Sigma }^{ab}`$ and $`R^{ab}`$ are considered as matrices. In the framework of the wave packet interpretation developed after Eq. (29), the function $`\mathrm{exp}(F_j(T))`$ plays the part of the space-time envelope of the wave packet associated with the jth neutrino mass eigenstate. The elements of the matrix $`Re(\mathrm{\Sigma }+i\frac{T}{2E_0}R)^1`$ constrain the extent of the wave packet envelope in space-time. As $`T`$ increases, the wave packet spreads because of the $`i\frac{T}{E_0}R`$ term. Thus the dispersion of the wave packet is due to the second order term in the expansion of the phase $`\varphi _j(𝐩_j)`$. Therefore, condition (37) or, equivalently, condition (38) means that dispersion has not yet begun in any direction, transversal or longitudinal. For that reason, the range of $`𝐋`$ values defined by $`Lp_0/\sigma _p^2`$ will be called the no-dispersion regime. Of course, this interpretation is not valid for $`\sigma _{x+}L`$, in which case the propagation time $`T`$ becomes indeterminate and dispersion loses its meaning.

Now that the origin of dispersion has been clarified, the term in Eq. (41) including $`R`$ can be neglected with respect to $`\mathrm{\Sigma }`$. Moreover, we approximate $`𝐯_j`$ by $`𝐯_0`$ in $`F_j(T)`$. In comparison with the oscillation formula (22) derived in section 3, this approximation will lead to the absence of the coherence-length term, since this term exclusively arises, when the dispersion is neglected, from the velocity difference $`𝐯_i𝐯_j`$. Dropping the index $`j`$, the wave packet envelope in Eq. (40) can then be written in the coordinate system diagonalizing $`\mathrm{\Sigma }`$ as

$$F(T)=\frac{(v_0^xTL^x)^2}{4\sigma _x^2}+\frac{(v_0^yTL^y)^2}{4\sigma _x^2}+\frac{(v_0^zTL^z)^2}{4\sigma _{x+}^2},$$
(42)
which shows that the conditions (36) assumed for Laplace’s method are required to obtain a non-negligible amplitude $`𝒜_j`$. If $`\sigma _{x+}L`$, condition (42) imposes the relation

$$𝐯_0T=𝐋+𝒪(\sigma _{x+}).$$

In that case, the phase of the amplitude (40) is given to $`𝒪(ϵ)`$ by

$$E_jT𝐩_j𝐋=E_0T𝐩_0𝐋+\frac{\delta m_j^2}{2p_0}\left(L+𝒪(\sigma _{x+})\right),$$
(43)
where the expansions (30) and (31) have been used. This result shows that the phase difference between the amplitude $`𝒜_i`$ and $`𝒜_j^{}`$ is equal to the standard oscillation phase present in Eq. (1).

What is the form of the amplitude (40) in the case of stationary boundary conditions ($`\sigma _{x+}L`$)? This question is important since such conditions were assumed by Ioannisian and Pilaftsis in their derivation of ‘plane wave oscillations’ . Let us first examine the wave packet envelope $`F(T)`$. In the stationary limit (see Eq. (34)), $`\sigma _{x+}`$ becomes larger than $`L^z`$ so that the $`z`$ component in Eq. (42) imposes that $`|v_0^z|T\sigma _{x+}`$. Putting this condition together with Eqs. (9) and (33), we obtain

$$|v_0^{x,y}|T\frac{|v_{\mathrm{SS}P,D}^{x,y}|}{|v_0^z|}\sigma _{x+}S_{\mathrm{SS}P,D}.$$

Thus the wave packet envelope (42) gives the following constraint in the stationary limit:

$$|L^{x,y}|S_{\mathrm{SS}P,D},$$
(44)
which means that, for propagation distances much larger than the size of the source and detector ($`LS_{\mathrm{SS}P,D}`$), the amplitude (40) is negligible unless the direction of observation $`𝐋`$ is nearly parallel to the average neutrino momentum $`𝐩_0`$. Thus this constraint is valid whether the stationary limit is taken or not.

Let us now examine the phase of the amplitude (40) in the stationary limit. The expansions (30) and (31) and the property $`\stackrel{~}{\rho }0`$ show that

$`E_j`$ $``$ $`E_0,`$
$`𝐯_0𝐩_j`$ $``$ $`𝐯_0𝐩_0{\displaystyle \frac{\delta m_j^2}{2E_0}}.`$
Using $`𝐋=𝐯_0L/|𝐯_0|`$ (with correction terms given by Eq. (44)), the phase of the amplitude is given to $`𝒪(ϵ)`$ by

$$E_jT𝐩_j𝐋=E_0T𝐩_0𝐋+\frac{\delta m_j^2}{2p_0}\left(L+𝒪(S_{\mathrm{SS}P,D})\right),$$
(45)
which leads again to the standard oscillation phase present in Eq. (1) since $`LS_{\mathrm{SS}P,D}`$. Eqs. (44) and (45) also show that apart from a factor $`exp(iE_0T)`$ which can be factorized from the sum over the mass eigenstates, the amplitude (40) is independent of time in the stationary limit.

We have thus demonstrated that the amplitude (40), derived under condition (38), is negligible in directions other than the neutrino average momentum whether the stationary limit is taken or not. Therefore the oscillation phase has the standard form given in Eq. (1) and no ‘plane wave oscillations’ can be observed, contrary to the claim made by Ioannisian and Pilaftsis . Also, the absence of the $`1/T`$ (or $`1/L`$) factor in Eq. (40), as noted by the same authors, is easily understood by noting that the absence of dispersion entails that the amplitude does not decrease with the distance.

#### 4.2.3 Probability

As in section 3, the last step towards the oscillation formula consists in computing the time average of the transition probability, which is a gaussian integral on $`T`$:

$$𝑑T𝒜_i𝒜_j^{}\mathrm{exp}\left(i\varphi _{ij}(T_0)\frac{(\delta m_i^2)^2+(\delta m_j^2)^2}{16\stackrel{~}{\sigma }_m^2E_0^2}\frac{(E_iE_j)^2}{4F^{\prime \prime }(T_0)}2F(T_0)\right),$$
(46)
where $`T_0`$ is the solution of $`F^{}(T_0)=0`$, while the phase $`\varphi _{ij}(T_0)`$ is given by

$$\varphi _{ij}(T_0)=(E_iE_j)T_0(𝐩_i𝐩_j)𝐋,$$
(47)
The function $`F(T)`$ can be rewritten as

$$F(T)=\frac{(\stackrel{~}{𝐯}_0T\stackrel{~}{𝐋})^2}{4\sigma _x^2},$$
(48)
where

$`\stackrel{~}{𝐯}_0`$ $`=`$ $`\sigma _x\sqrt{\mathrm{\Sigma }^1}𝐯_0,`$
$`\stackrel{~}{𝐋}`$ $`=`$ $`\sigma _x\sqrt{\mathrm{\Sigma }^1}𝐋,`$
with $`\mathrm{\Sigma }^1=diag(\sigma _x^2,\sigma _x^2,\sigma _{x+}^2)`$. It is then easy to compute

$`T_0`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{𝐯}_0\stackrel{~}{𝐋}}{\stackrel{~}{𝐯}_0^2}},`$ (49)
$`F(T_0)`$ $`=`$ $`{\displaystyle \frac{(\stackrel{~}{𝐯}_0\times \stackrel{~}{𝐋})^2}{4\sigma _x^2\stackrel{~}{𝐯}_0^2}},`$ (50)
$`F^{\prime \prime }(T_0)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{𝐯}_0^2}{2\sigma _x^2}}.`$ (51)
We now show that results (43)-(45), obtained at the level of the amplitude, can be reproduced by a careful analysis of Eq. (46). This is not completely trivial as the integration variable $`T`$ becomes indeterminate in the stationary limit.

The insertion of expression (50) in Eq. (46) shows that $`F(T_0)`$ plays the role of a directional constraint. More precisely, $`𝑑T𝒜_i𝒜_j^{}`$ is non-negligible if $`\stackrel{~}{𝐋}`$ is nearly parallel to $`\stackrel{~}{𝐯}_0`$ or, for $`\sigma _{x+}L`$, if $`𝐋`$ is nearly parallel to $`𝐯_0`$:

$$𝐋=\frac{𝐯_0}{|𝐯_0|}L+𝒪(\sigma _{x+}).$$
(52)
With the substitution (52), the value of $`T_0`$ given by Eq. (49) becomes $`T_0=L/|𝐯_0|+𝒪(\sigma _{x+})`$, so that the phase (47) becomes

$$\varphi _{ij}(T_0)=\frac{\delta m_{ij}^2}{2p_0}(L+𝒪(\sigma _{x+})),$$
(53)
which is equal to the standard oscillation phase present in Eq. (1) and equivalent to result (43).

If $`\sigma _{x+}L`$, $`\stackrel{~}{𝐯}_0`$ and $`\stackrel{~}{L}^z`$ tend to zero so that $`F(T_0)`$ should directly be studied as a function of $`𝐯_0`$ and $`𝐋`$. General conclusions for arbitrary $`\sigma _{x\pm }`$ can be drawn from the study of the quadratic form in $`(L^x,L^y,L^z)`$ associated with $`F(T_0)=1`$. This analysis shows that there is an eigenvalue $`s_3=0`$ corresponding to an eigenvector along $`𝐯_0`$. The two other eigenvalues $`s_{1,2}`$ are positive (with $`s_2s_1`$) , so that the surface $`F(T_0)=1`$ in $`(L^x,L^y,L^z)`$-space is a cylinder of elliptical section with an axis along $`𝐯_0`$. This geometrical picture can be interpreted as imposing that the components of $`𝐋`$ orthogonal to $`𝐯_0`$ should be smaller than $`\sqrt{1/s_2}`$, whereas there is no constraint at all on the component of $`𝐋`$ along $`𝐯_0`$. In the stationary limit, the lengthy expressions of the non-zero eigenvalues become

$`s_1`$ $``$ $`{\displaystyle \frac{1}{4\sigma _x^2}},`$
$`s_2`$ $``$ $`{\displaystyle \frac{1}{4\sigma _x^2}}{\displaystyle \frac{\sigma _x^2(v_0^z)^2}{\sigma _{x+}^2(v_0^x)^2+\sigma _{x+}^2(v_0^y)^2+\sigma _x^2(v_0^z)^2}}.`$ (54)
The properties of the stationary limit, namely Eqs. (9), (33) and (34), lead to the bound $`\sqrt{1/s_2}S_{\mathrm{SS}P,D}`$. In the stationary limit, the components $`L^{x,y}`$ (which are then orthogonal to $`𝐯_0`$) should thus be smaller than $`S_{\mathrm{SS}P,D}`$. Thus the time-averaged probability is always negligible in directions other than along the average neutrino momentum $`𝐩_0`$. Note that this condition reproduces the constraint (44) derived in the study of the stationary limit of the amplitude (40). The phase $`\varphi _{ij}(T_0)`$ can then be evaluated as in Eq. (45), leading to the oscillation phase

$$\varphi _{ij}(T_0)=\frac{\delta m_{ij}^2}{2p_0}(L+𝒪(S_{\mathrm{SS}P,D})),$$
(55)
which is again equal to the standard oscillation phase present in Eq. (1) and equivalent to result (45).

Putting together the results (46), (50), (51) and (53) (or (55)), the flavor-mixing transition probability (5) for a propagation distance satisfying condition (38) can be written as

$`𝒫_{\alpha \beta }(𝐋)\mathrm{exp}\left({\displaystyle \frac{(\stackrel{~}{𝐯}_0\times \stackrel{~}{𝐋})^2}{2\sigma _x^2\stackrel{~}{𝐯}_0^2}}\right){\displaystyle \underset{i,j}{}}V_{i\alpha }V_{i\beta }^{}V_{j\alpha }^{}V_{j\beta }`$
$`\times \mathrm{exp}\left(2\pi i{\displaystyle \frac{L}{L_{ij}^{osc}}}2\pi ^2\left({\displaystyle \frac{\stackrel{~}{\rho }\stackrel{~}{\sigma }_{xeff}}{L_{ij}^{osc}}}\right)^2{\displaystyle \frac{(\delta m_i^2)^2+(\delta m_j^2)^2}{16\stackrel{~}{\sigma }_m^2E_0^2}}\right),`$ (56)
where $`\stackrel{~}{\sigma }_{xeff}=\sigma _x|𝐯_0|/|\stackrel{~}{𝐯}_0|`$ and $`L_{ij}^{osc}`$ is given by Eq. (23).

This oscillation formula is very similar to Eq. (22). There are three main differences, two of which have already been mentioned. First there is no coherence-length term in Eq. (56), which is due to the neglect of the terms beyond $`𝒪(ϵ^0)`$ in the evaluation of $`F(T_0)`$ (see Eq. (42)). Secondly, the geometrical decrease in $`1/L^2`$ is lacking in Eq. (56), which is explained by the fact that the dispersion of the neutrino ‘wave packet’ is not yet significant. Finally, the directional constraint present in Eq. (22) restricts the neutrino propagation to a cone of angle $`arcsin(\sigma _p/p_0)`$, whereas the directional constraint in Eq. (56) confines the propagation to a cylinder of radius $`r`$ (with $`\sigma _xrS_{\mathrm{SS}P,D}`$). This different behavior is also a result of the absence of dispersion for $`Lp_0/\sigma _p^2`$.

Before closing this section, it is interesting to understand why a directional constraint is missing in Ref.  for $`Lp_0/\sigma _p^2`$, as this fact explains the ‘plane wave oscillations’ result. At the end of their computations, the authors of Ref.  obtain an amplitude $`𝒜_j`$ whose dominant term depends on $`exp(ip_j|𝐋^{}|)`$, where $`|𝐋^{}|=\sqrt{𝐋_{}^{}{}_{}{}^{2}}`$ is the complex ‘norm’ of a complex vector $`𝐋^{}=𝐋2i\sigma _x^2𝐩_0`$ (we have translated their results in our notation through the correspondence $`q_jp_j`$, $`\stackrel{}{k}𝐩_0`$, $`\delta l^24\sigma _x^2`$, $`\stackrel{}{L}𝐋^{}`$, $`\stackrel{}{l}𝐋`$). For $`Lp_0/\sigma _p^2`$, the quantity $`|𝐋^{}|`$ can be expanded and the argument of the exponential reads to second order

$$ip_j|𝐋2i\sigma _x^2𝐩_0|2\sigma _x^2p_0p_j+ip_j\frac{𝐩_0𝐋}{p_0}\frac{p_j}{4\sigma _x^2p_0}\left(L^2\frac{(𝐩_0𝐋)^2}{p_0^2}\right).$$

The two last terms were neglected in Ref.  though $`L\sigma _x`$. They lead to the following directional constraint:

$$\mathrm{exp}\left(\frac{(𝐯_0\times 𝐋)^2}{4\sigma _x^2v_0^2}\right),$$
(57)
where the factor 4 becomes a factor 2 when the amplitude is squared. Therefore the result of Ref.  also includes for $`Lp_0/\sigma _p^2`$ a directional constraint forbidding plane wave oscillations.

Note that the stationary limit assumption $`|𝐯_{\mathrm{SS}P,D}|/\sigma _{e\mathrm{SS}P,D}=0`$ leads to $`s_1=s_2`$ (see Eq. (54)), so that the constraint (57) becomes equal to the one present in Eq. (56). Thus the condition $`|𝐯_{\mathrm{SS}P,D}|/\sigma _{e\mathrm{SS}P,D}=0`$ seems generic for models with stationary boundary conditions, since it was also applied in the case of the Grimus-Stockinger model (see discussion after Eq. (9) and at the end of section 3).

In conclusion, the analysis of the amplitude (40) (as well as of the transition probability (56)) derived in the external wave packet model under condition (38) shows that ‘plane wave oscillations’ do not exist. This negative result is confirmed by a reexamination of the formula derived in Ref. . Besides providing a rigorous method to compute the transition probability (56), the external wave packet model has the advantage of associating a clear physical picture to the different stages of computations. For example, condition (38) can be interpreted as meaning that the wave packet associated with the neutrino has not yet begun to spread in any direction. Finally, the oscillation formula (56) reduces to the standard formula (1) if the observability conditions $`|𝐯_0\times 𝐋|v_0\sigma _x`$ and $`L_{ij}^{osc}\sigma _x`$ are satisfied.

### 4.3 Transversal- and longitudinal-dispersion regimes

Let us now assume that condition (38) is not satisfied, i.e. $`Lp_0/\sigma _p^2`$. In that case, Laplace’s method cannot be used to integrate on all three components $`p^{x,y,z}`$ in the amplitude (29) since dispersion becomes significant. However the spreading of the neutrino ‘wave packet’ is not identical in all directions. More specifically, the onset of dispersion in the direction $`𝐩_0`$ can be delayed by two factors. First, the matrix element $`R^{ab}`$ present in Eq. (35) leads to a relativistic contraction (of $`1𝐯_0^2`$) in the direction $`𝐩_0`$ of the dispersion of the neutrino ‘wave packet’ (see Eq. (41)). Secondly, the momentum width along $`𝐩_0`$ is given for $`\sigma _{e\mathrm{SS}P,D}\sigma _x`$ (i.e. in the stationary limit (34)) by a vanishing $`\sigma _{p+}`$. Thus Laplace’s method will be valid for a longer time $`T`$ in the direction $`𝐩_0`$ than in directions transverse to this vector. For this reason, the choice of the integration method in the direction $`𝐩_0`$ will be postponed for a short while, while the method of stationary phase will be preferred for momentum integrations in directions transverse to $`𝐩_0`$.

Let the $`z`$ axis be along $`𝐋`$, i.e. $`𝐋=L𝐞_z`$. As in section 4.2, the examination of the amplitude (29) shows that the quick variation of the phase averages the amplitude to zero unless $`𝐩_0`$ and $`𝐋`$ are nearly parallel. The method of stationary phase will thus be applied in directions $`p^{x,y}`$, the stationary points of which are given by $`p^x=p^y=0`$. The result of the method of stationary phase for the transverse momenta in (29) can be written as follows:

$$𝒜_j\frac{g(𝐥)}{Ti\mu }𝑑pe^{i\varphi _j(p)f_j(p)},$$
(58)
where $`pp^z`$ and $`\mu =E_0/2\sigma _p^2`$. The geometrical constraint $`g(𝐥)`$ is given by Eq. (14), while the momentum width $`\sigma _p`$ is defined by Eq. (15). The phase $`\varphi _j(p)`$ is given by

$$\varphi _j(p)=\sqrt{p^2+m_j^2}TpL.$$
(59)
The original overlap function is partly included in $`g(𝐥)`$ and partly in

$$f_j(p)=f_{j\mathrm{SS}P}(p)+f_{j\mathrm{SS}D}(p),$$
(60)
with

$$f_{j\mathrm{SS}P}(p)=\frac{(pp_0)^2}{4\sigma _{p\mathrm{SS}P}^2}+\frac{\left(\sqrt{p^2+m_j^2}E_0(pp_0)v_{\mathrm{SS}P}\right)^2}{4\sigma _{e\mathrm{SS}P}^2},$$

where $`v_{\mathrm{SS}P}=v_{\mathrm{SS}P}^z`$ and $`p_0=p_0^z`$. The energy $`E_0`$ has been redefined so as to absorb a factor $`p_0^xv_{\mathrm{SS}P}^x+p_0^yv_{\mathrm{SS}P}^y`$. The definition of $`f_{j\mathrm{SS}D}(p)`$ is similar. As $`TE_0/\sigma _{peff}^2`$, the term $`1/(Ti\mu )`$ can be approximated by $`1/T`$ so as to provide a prefactor $`1/L^2`$ in the transition probability. This expected geometrical decrease is seen to originate in the transverse dispersion of the wave packet corresponding to the oscillating particle. For future use, we define a reference mass $`m_0`$ by $`m_0^2=E_0^2p_0^2`$ and a velocity $`v_0`$ by $`v_0=p_0/E_0`$.

As regards the longitudinal momentum integral, it is tempting to proceed as in section 3, i.e. to compute first the time average of the transition probability before integrating on $`p`$. However, the prefactor $`1/(Ti\mu )`$ present in the amplitude (58) gives in that case a delta function look-alike of width $`\mu ^1\sigma _p^2/E_0`$, introducing an additional momentum uncertainty which is larger than the mass difference $`\delta m^2/E`$ since $`L_{ij}^{osc}p_0/\sigma _p^2`$. For this reason, it is preferable to avoid this shortcut (note that it yields the same final answer as given by the following method). Moreover, it is interesting for the physical interpretation to postpone the time average, so as to obtain the explicit dependence of the amplitude on time and distance. As in the case of transverse momenta integrals, the choice of the method to perform the longitudinal momentum integration (58) is done by comparing the expansions of the phase and of the overlap function around the value $`p_j`$ for which $`f_j(p)`$ is extremal. One obtains to $`𝒪(ϵ)`$

$$p_j=p_0+(\rho 1)\frac{\delta m_j^2}{2p_0},$$

where the dimensionless number $`\rho `$ is defined by Eq. (16). The associated energy $`E_j=\sqrt{p_j^2+m_j^2}`$ and velocity $`v_j=p_j/E_j`$ are given to $`𝒪(ϵ)`$ by

$`E_j`$ $`=`$ $`E_0+\rho {\displaystyle \frac{\delta m_j^2}{2E_0}},`$
$`v_j`$ $`=`$ $`v_0+\left(\rho (1v_0^2)1\right){\displaystyle \frac{\delta m_j^2}{2p_0E_0}}.`$ (61)
The expansions of $`f_j(p)`$ and $`\varphi _j(p)`$ are given by

$`f_j(p)`$ $``$ $`{\displaystyle \frac{(pp_j)^2}{4\sigma _{peff}^2}},`$
$`\varphi _j(p)`$ $``$ $`\varphi _j(p_j)+(v_jTL)(pp_j)+{\displaystyle \frac{m_j^2T}{2E_0^3}}(pp_j)^2,`$ (62)
where $`\sigma _{peff}`$ is given by Eq. (17) and can again be interpreted as the longitudinal width of the overlap function. Laplace’s method will be used if the phase (62) varies slowly over the width $`\sigma _{peff}`$, i.e. if the two following conditions are satisfied

$`|v_jTL|2\sigma _{peff}`$ $``$ $`1,`$ (63)
$`{\displaystyle \frac{m_j^2T}{2E_0^3}}4\sigma _{peff}^2`$ $``$ $`1.`$ (64)
As in section 4.2, the first order constraint (63) is included in the result of Laplace’s method. Thus the criterion allowing to choose between Laplace and stationary phase methods is given by Eq. (64). In other words, it is better to use Laplace’s method if $`T`$ is smaller than a dispersion time $`T_j^{disp}`$ defined by

$$T_j^{disp}=\frac{E_0^3}{2m_j^2\sigma _{peff}^2}.$$
(65)
The term ‘dispersion time’ is justified by the fact that it is the time at which the longitudinal dispersion of the amplitude becomes important, more precisely twice the initial size. A dispersion length $`L_j^{disp}`$ can be defined by $`L_j^{disp}=v_0T_j^{disp}`$. The distance range $`p_0/\sigma _{peff}^2LL_j^{disp}`$ will be called the tranverse-dispersion regime. For $`LL_j^{disp}`$, the stationary phase method is more accurate: this distance range will be called the longitudinal-dispersion regime. Various estimates of the dispersion length are given in , showing that the concept of dispersion length is relevant to nonrelativistic particles as well as to supernova neutrinos, and possibly to solar neutrinos.

In the transverse-dispersion regime, the evaluation of the amplitude (58) as a gaussian integral around $`p_j`$ gives

$$𝒜_j\frac{g(𝐥)\sigma _{peff}}{T\sqrt{1+iT/T_j^{disp}}}\mathrm{exp}\left(iE_jT+ip_jL\left(\frac{\delta m_j^2}{4\sigma _mE_0}\right)^2\frac{1}{1+iT/T_j^{disp}}\frac{(v_jTL)^2}{4\sigma _{xeff}^2}\right).$$
(66)
The width $`\sigma _m`$ is defined by Eq. (20). The amplitude (66) behaves as a wave packet of group velocity $`v_j`$ and space-time extent $`(1+(T/T_j^{disp})^2)^{\frac{1}{2}}\sigma _{xeff}`$. If the longitudinal dispersion is neglected ($`T_j^{disp}=\mathrm{}`$), the amplitude (66) is similar to Eq. (18) of Ref. .

The following step consists in computing the time averaged transition probability (5). If dispersion is taken into account ($`T_j^{disp}<\mathrm{}`$), a tedious gaussian integration by Laplace’s method (see ) yields the oscillation formula (22). If the dispersion is neglected ($`T_j^{disp}=\mathrm{}`$), a much shorter computation leads to the same oscillation formula (22), except that the following substitution has to be made:

$$L_{ij}^{coh}\frac{L_{ij}^{coh}}{|\rho (1v_0^2)1|}(FALSE).$$
(67)
The incorrect multiplying factor has its origin in Eq. (61): in the limit $`T_j^{disp}=\mathrm{}`$, the coherence length term arises only from the difference between the group velocities $`v_i`$ and $`v_j`$. However the factor $`|\rho (1v_0^2)1|`$ tends to 1 in the relativistic limit, so that the substitution (67) becomes trivial. This observation explains why our result (22) coincides with Eq. (26) of Ref. , as the authors of this article, while neglecting dispersion, consider only relativistic neutrinos. Note however that even relativistic neutrinos spread at large distances so that a calculation neglecting dispersion such as in Ref.  is only valid for $`LL_j^{disp}`$.

At sufficiently large distance, dispersion becomes significant and all neutrinos propagating freely enter into the longitudinal-dispersion regime. In this regime, we have argued that the integral (58) should be evaluated with the method of stationary phase. The stationary point of the phase $`\varphi _j(p)`$ is given by

$$p_{cl,j}=m_j\frac{v_{cl}}{\sqrt{1v_{cl}^2}},$$

where $`v_{cl}=L/T`$. It can be interpreted as the classical momentum of a particle of mass $`m_j`$, travelling at the classical velocity $`v_{cl}`$. The evaluation of the amplitude (58) as a gaussian integral around $`p_{cl,j}`$ gives

$$𝒜_j\frac{g(𝐥)\sigma _{peff}}{T\sqrt{1+iT/T_j^{disp}}}\mathrm{exp}\left(im_j\sqrt{T^2L^2}f_j(p_{cl,j})+\sigma _{peff}^2\frac{\left(f_j^{}(p_{cl,j})\right)^2}{1+iT/T_j^{disp}}\right),$$
(68)
where $`f_j(p)`$ is defined by Eq. (60), $`f_j^{}(p)`$ refers to its derivative. The wave packet interpretation of the amplitude (68) is not obvious but the shape of the associated wave packet can be studied by an expansion around the maximum of the amplitude. Once more, the time average of the transition probability (5) is computed by Laplace’s method (see Ref. ), yielding again the oscillation formula (22).

It is striking that the two different methods of approximation used in this section, a priori valid in different regimes, give the same oscillation formula, which also coincides with the result obtained with the approximation method used in section 3, Eq. (22). The dispersion length does not play any special role in the final result. Thus, each method is accurate enough to be extended to the whole range of distances. However, the physical interpretation of the terms appearing in the oscillation formula (22) depends on the relative values of $`L`$ and $`L_{i,j}^{disp}`$. From the detailed analysis done in , these terms can be interpreted as follows:

* In the transverse-dispersion regime, the coherence-length term mainly arises from the progressive spatial separation of ‘wave packets’, due to their different group velocities. The localization term mainly arises from the initial spread of the ‘wave packet’, which must be smaller than the oscillation length. This is in agreement with the traditional explanation of these two terms. Note that $`L_i^{disp}L_j^{disp}`$ if $`m_im_j`$. In that case, one can check that decoherence occurs for $`L<L_i^{disp}`$.
* In the longitudinal-dispersion regime, the coherence-length term mainly arises from the spreading of the neutrino ‘wave packet’: the interference term is averaged to zero by the time integral as the spread of the wave packet becomes large. The localization term mainly arises from the separation of the ‘wave packets’: if they are spatially separate at the onset of the longitudinal-dispersion regime, they will never again overlap because their separation increases as fast as their dispersion. This is a new explanation for the origin of these two terms.

Thus the origins of the coherence-length and localization terms in the transverse-dispersion regime become respectively the origins of the localization and coherence-length terms in the longitudinal-dispersion regime. The situation can be summarized in the following diagram, where ‘wp’ is an abbreviation for ‘wave packet’:

|  | transverse-dispersion regime | longitudinal-dispersion regime |
| --- | --- | --- |
| Coherence length | decreasing overlap of wp | increasing dispersion of wp |
| Localization | initial spread of wp | constant overlap of wp |

TABLE I. Origins of the coherence-length and localization terms.

It is now possible to examine Shtanov’s analysis in configuration space , since the explicit $`L,T`$ dependence of the amplitude $`𝒜_j`$ has been computed for the whole range of distances (see Eqs. (40), (66) and (68)). This author works directly with the neutrino propagator in configuration space, and considers the space-time variables $`x,x^{}`$ (see FIG. 1) as macroscopic variables. His model will be considered in a scalar version in order to make the comparison with our results easier. In configuration space, the propagator for a scalar neutrino propagating from $`x`$ to $`x^{}`$ is asymptotically ($`m_j\sqrt{(x^{}x)^2}1`$) given by

$$𝒜_j\frac{\sqrt{m_j}}{\left((x^{}x)^2\right)^{3/4}}e^{im_j\sqrt{(x^{}x)^2}},$$
(69)
where $`(x^{}x)^2`$ is the Lorentz interval. Shtanov then computes the convolution of this propagator with a source. In the case of a monochromatic source, he obtains the standard oscillation formula. However, a convolution with a strongly localized source ($`\sigma _{x,t}1/E`$) leads to an amplitude that keeps its dependence on the mass prefactor $`\sqrt{m_j}`$. In that case, the transition probability is not equivalent to the oscillation formula (22), at least if the neutrino masses are not nearly degenerate. We proceed to show that Shtanov’s result is incorrect. Note first that the amplitude (69) is in correspondence with our amplitude (68), computed with the stationary phase method, since the prefactor in Eq. (68) can be rewritten, for $`TT_j^{disp}`$, as

$$\frac{1}{T}\frac{\sigma _{peff}}{\sqrt{1+iT/T_j^{disp}}}\frac{\sqrt{m_j}}{(T^2L^2)^{3/4}},$$

where the definition (65) of the dispersion time has been used. This prefactor coincides with the prefactor in Eq. (69). However, the subsequent time average of the transition probability completely cancels this dependence on the mass, yielding Eq. (22). This can be seen by expanding the argument of the exponential in Eq. (68) around the average propagation time. The width with respect to $`T`$ of the amplitude $`𝒜_j`$ is found to be equal to $`\frac{\sigma _{peff}T}{m_jv_0\gamma _{cl}^3}`$ (where $`\gamma _{cl}`$ is the Lorentz factor associated with the velocity $`v_{cl}=L/T`$), thereby providing a $`m_j`$-dependent factor that cancels the $`m_j`$ prefactor in the gaussian integration.

Shtanov does not obtain such a result, since he does not average the probability on time. Besides, the mass prefactors would remain even if a time average were performed on the probability: no ‘wave packet envelope’ appears indeed in Shtanov’s amplitude, as this author does not try to compute the explicit convolution of the propagator (69) with a wave packet source of arbitrary width. For that reason, he wrongly concludes that the mass prefactor remains if the source is strongly localized. Shtanov also derives another oscillation formula for mass eigenstates such that $`m_2\sqrt{(x^{}x)^2}1m_1\sqrt{(x^{}x)^2}`$. However, it can be checked that decoherence occurs in that case. In conclusion, Shtanov’s computations in configuration space do not lead to new oscillation formulas.

The last case to be considered is the oscillation formula obtained by Blasone, Vitiello and collaborators (hereafter abbreviated as BV). These authors have attempted to define a Fock space of weak eigenstates and to derive a nonperturbative oscillation formula . They define flavor creation and annihilation operators, satisfying canonical (anti)commutation relations, by means of Bogoliubov transformations. Apart from the speculative nature of the undertaking, the drawbacks of the approach are the dependence on time, not on space, of the oscillation formula, as well as the lack of observability conditions. The latter problem is important since it determines whether the new features of the BV oscillation formula are observable in practice. Since the BV oscillation formula tends to the oscillation formula (22) in the relativistic limit or if the mass eigenstates are nearly degenerate (with the coherence-length and localization terms removed), we can focus on the case of a nonrelativistic oscillating particle having very distinct mass eigenstates. In that case, $`p\delta m^2/2E`$, so that either $`\sigma _p\delta m^2/2E`$ or $`p\sigma _p`$. Under these conditions, the oscillation formula (22) shows that oscillations vanish, either because of the localization term or because of the coherence-length term (though Eq. (22) is strictly speaking not valid in the case considered, localization and coherence-length conditions are generic features of quantum-mechanical models of oscillations ). Once the oscillation terms have been averaged to zero both in the BV formula and in Eq. (22), these two oscillation formulas do not differ anymore. Therefore, the BV formalism does not seem to be relevant to the phenomenology of oscillations on macroscopic distances. This observation does not detract from the theoretical worth of that approach.

We shall not treat here the charged lepton oscillations obtained by Srivastava, Widom and Sassaroli . The reader is referred to  for a discussion of this subject.

## 5 Conclusion

In this article, we have shown that all correct results obtained by applying quantum field theory to neutrino oscillations far from the source are incorporated in the oscillation formula (22). Closer to the source (i.e. for $`Lp_0/\sigma _p^2`$), this formula is modified in only one respect without repercussions on the oscillation length: the global spatial distribution of the detection probability has the form given in Eq. (56), with the consequence that the neutrino propagation is restricted within a cylinder instead of a cone. The quantum field oscillation formulas (22) and (56) are equivalent to the standard result (1) if some observability conditions are verified. The new oscillations formulas proposed by Ioannisian and Pilaftsis  and by Shtanov  have been disproved, while the Blasone-Vitiello formula  has been shown to be phenomenologically equivalent to the standard result. Besides, arguments against oscillations of charged leptons, as proposed by Srivastava, Widom and Sassaroli , can be found in Ref. .

We have insisted on the fact that the flavor-mixing transition amplitude can be interpreted in terms of wave packets, so that oscillations can be seen, like in the intermediate wave packet model, as the result of an interference between propagating wave packets. Moreover, the equivalence (11) shows that this physical picture still holds in the case of stationary boundary conditions, provided that an incoherent sum over the energy is performed, contrary to what was claimed in Ref. . As a byproduct of our analysis, the propagation of a particle has been shown to go through three regimes, separated by two thresholds which are marked by the onset of first the transversal dispersion and then the longitudinal dispersion of the associated wave packet. This picture has been found to be useful for the understanding of the terms present in the oscillation formula as well as for the computation of the correct spatial distribution of the probability. We believe that the wave packet interpretation is the only consistent picture of the oscillation process. In particular, it shows the irrelevance of discussions bearing on the neutrino energy-momentum, since it has not a unique value for a given mass eigenstate, but rather a spread described by the overlap function $`\psi (\sqrt{𝐩^2+m_j^2},𝐩)`$.

One could wonder what advantage the external wave packet model has over the intermediate wave packet model, apart from being the only consistent way to derive the oscillation formula. Actually, the quantum field formula (22) has two important new features with respect to the quantum mechanical result. First, the energy uncertainty is not put in by hand, but is defined in terms of the 3-momentum widths and velocities of the in- and outgoing wave packets (see Eq. (7)), with the result that the stationary limit is well-defined and that it is possible to analyze the dependence of the coherence length on the accuracy of energy measurements. Secondly, the oscillation formula is valid as well for neutrinos as for K and B mesons, putting it on a much firmer phenomenological basis. This might be useful in the light of the numerous nonstandard oscillation formulas existing in the literature (see Ref.  for a review).

The stationary limit has been considered in our computations to show the great generality of the external wave packet model. The stationary assumption is not problematic for the derivation of the oscillation formula. However, this limiting case is rather unphysical. The computations of section 4 have indeed shown that boundary conditions can be considered as stationary when $`\sigma _{xeff}L`$. This constraint, which can be rewritten either as $`\sigma _{e\mathrm{SS}P,D}1/L`$ or as $`T_{\mathrm{SS}P,D}^{overlap}T`$ (verified either at the source or at the detector), is extremely stringent. In the example of atmospheric neutrinos, the process can be considered as stationary if $`\sigma _{e\mathrm{SS}P,D}10^{19}`$MeV or equivalently if $`T_{\mathrm{SS}P,D}^{overlap}10^3`$s. These conditions are far from being realized at the microscopic level in any known neutrino source or detector. Although most neutrino sources are stationary from a macroscopic point of view, there is no reason to think that individual particles in the source and detector remain unperturbed in coherent states over macroscopic time scales: the Sun is certainly not stationary at the atomic scale, and neither is a detection process where charged leptons are observed with a finite energy and time spread .

A subject deserving further investigation in the quantum field framework is the influence of the decay width of the source of the neutrino. The case of decay at rest has been studied in detail with the Wigner-Weisskopf approximation in Ref. , while the treatment of decay in flight  is still unsatisfactory. Unstable neutrinos have been studied in Ref. .

We would be sorry to leave the reader with the impression that the quantum field approach to oscillations is complicated. Actually, this is a consequence of the generality of the external wave packet model. For pedagogical purposes, it is interesting to resort to a simplified quantum field model, where the source and detector are perfectly located in space and in- and outgoing states are stationary .

## Acknowledgments

I thank Jeanne De Jaegher for helpful comments on the manuscript.