1 Introduction

The idea of neutrino oscillations as a mechanism to solve the solar-neutrino puzzle was firstly proposed by Pontecorvo and collaborators [1,2,3,4] and it was later confirmed by a plethora of experiments (see e.g. [5,6,7,8]).

Although many features of neutrino mixing and oscillations are now well-understood [9,10,11], there is no agreement on their correct ultimate description within quantum field theory (QFT). Various ideas were proposed in the last three decades, as external wavepackets [12, 13], weak-process states [14] and the flavor Fock-space approach [15,16,17,18,19,20,21]. The latter is based on the discovery [15] that the flavor and the mass representations of the equal-time anticommutation relations of neutrino fields, are unitarily inequivalent [22,23,24,25]. Therefore, the Hilbert space where flavor fields are defined is explicitly built in and the oscillation probability is computed by taking the expectation value of lepton currents/charges on the one-particle neutrino states at a reference time. Such modified formula differs from the classic Pontecorvo result in two respects [26]: (i) apart from the usual oscillation term which depends on the difference of neutrino energies/frequencies, the oscillation formula of Ref. [26] shows up a fast-oscillation term which depends on the sum of the frequencies; (ii) in the formula of Ref. [26], there are energy dependent oscillation amplitudes which are the coefficients of a Bogoliubov transformation [15].

In this paper we introduce a different approach, in a close analogy to what is done in the study of unstable particles [27, 28]. In fact, we employ the interaction (Dirac) picture, where the interaction Lagrangian in the Dyson series only contains the mixing term between different flavor-fields. For simplicity, we limit our calculation to the case of two flavors. Then we compute amplitudes for the various decay channels at the first order, which describe both flavor changing and survival processes. Three examples are here analyzed: a quantum mechanical (QFT in 0+1D) toy model, a scalar field model and a fermion (“neutrino”) model in 3+1D. Remarkably, we find the that the fermion flavor-transition formula non-trivially agrees, within the approximation adopted, with the non-perturbative formula of the flavor-Fock space approach. Let us remark that the comparison is not possible in the boson case, where the flavor charge expectation value is not positive-definite and thus it cannot be interpreted as a probability [29, 30]. In this respect, the present work represents also a viable approach to compute the oscillation probability in this tricky situation.

The paper is organized as follows: in Sect. 2 we present general considerations on field mixing and the interaction picture approach. In Sect. 3 we study the 0+1D toy model, while in Sects. 4 and 5 we extend our consideration to 3+1D scalar and fermion models, respectively. Finally, we present discussion and conclusions in Sect. 6. For reader’s convenience, in Appendix A we briefly review the non-perturbative flavor Fock space approach.

2 General considerations

The charged-current lepton sector of weak interaction is described (in the case of two-flavors) by the Lagrangian

(1)

with

$$\begin{aligned} \mathcal {L}_{mix}= & {} -m_{e \mu } \left( \overline{\nu }_e \nu _\mu +\overline{\nu }_\mu \nu _e\right) \,, \end{aligned}$$
(2)
$$\begin{aligned} \mathcal {L}_{wint}= & {} -\frac{g}{2\sqrt{2}}\sum _{\sigma =e,\mu }\left[ W_{\mu }^{+}\, \overline{\nu }_\sigma \,\gamma ^{\mu }\,(1-\gamma ^{5})\,l_\sigma + h.c. \right] \end{aligned}$$
(3)

The neutrino kinetic term (including \(\mathcal {L}_{mix}\)) can be diagonalized by the mixing transformation [31, 32]

$$\begin{aligned} \nu _\sigma \ = \ \sum _{j=1,2} U_{\sigma j}^* \nu _j \,, \end{aligned}$$
(4)

U is the mixing matrix. In the two flavor case, here analyzed

$$\begin{aligned} U \ = \ \begin{pmatrix} \cos \theta &{} \sin \theta \\ -\sin \theta &{} \cos \theta \end{pmatrix} \,, \end{aligned}$$
(5)

with \(\tan 2 \theta = 2 m_{e \mu }/(m_\mu -m_e)\).

If one employs the interaction picture to compute transition amplitudes, \(\mathcal {L}\) must be decomposed into a free and an interaction part. A possible choice is

$$\begin{aligned} \mathcal {L} \ = \ \mathcal {L}_0^m+\mathcal {L}_{int}^m \,, \end{aligned}$$
(6)

with

$$\begin{aligned} {{\mathcal {L}}}_0^m= & {} \sum _j \, \overline{\nu }_j \left( i \gamma _\mu \partial ^\mu - m_j \right) \nu _j \, \nonumber \\ {}{} & {} \quad + \, \sum _\sigma \, {\overline{l}} \left( i \gamma _\mu \partial ^\mu - {\tilde{m}}_{\sigma } \right) l \,, \end{aligned}$$
(7)
$$\begin{aligned} {{\mathcal {L}}}_{int}^m= & {} -\frac{g}{2\sqrt{2}} \sum _{\sigma ,j} \,\left[ W_{\mu }^{+}\, \overline{\nu }_{j} \, U^*_{j \sigma }\,\gamma ^{\mu }\,(1-\gamma ^{5})\, l_\sigma + h.c. \right] \,. \end{aligned}$$
(8)

In such a case the effect of mixing is incorporated in the weak-interaction vertex. Following this approach, one is led to calculate transition amplitudes in which neutrinos appear only as internal lines [12, 13, 33].Footnote 1

However, in charged current weak interaction processes, neutrinos are produced with a definite flavor. Therefore, another reasonable possibility is to take the following split

$$\begin{aligned} \mathcal {L} \ = \ \mathcal {L}_{0}+\mathcal {L}^{g}_{int} \,, \end{aligned}$$
(9)

with

(10)
$$\begin{aligned} \mathcal {L}^g_{int}= & {} \mathcal {L}_{mix}+\mathcal {L}_{wint} \,. \end{aligned}$$
(11)

In this approach, \(\mathcal {L}_{wint}\) is diagonal in the asymptotic fields appearing in Eq. (10). Thus, in order to describe neutrino oscillations, we can safely disregard \(\mathcal {L}_{wint}\) (zeroth-order in g), so that the charged-lepton part also decouples. In other words, we can treat the mixing term as an interaction, and we can compute the transition amplitudes among different flavors by means of the usual Dyson formula for the time evolution operator

$$\begin{aligned} U(t_i,t_f) \ {}= & {} \ \mathcal {T} \exp \left[ i \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x \,:\mathcal {L}_{int}(x): \right] \nonumber \\= & {} \ \mathcal {T}\exp \left[ -i \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x \,:\mathcal {H}_{int}(x): \right] \,, \end{aligned}$$
(12)

where \(\mathcal {L}_{int} \equiv \mathcal {L}^{g=0}_{int}=-m_{e \mu } \left( \overline{\nu }_e \nu _\mu +\overline{\nu }_\mu \nu _e\right) \), \(\mathcal {H}_{int}(x)=-\mathcal {L}_{int}(x)\) is the interaction Hamiltonian density and \(\mathcal {T}\) is the chronological product. In the following we will only need the expression of the operator up to the second order

$$\begin{aligned} U(t_i,t_f) \ {}= & {} \ 1-i\int _{t_{i}}^{t_{f}}\!\!\textrm{d}t_{1} \,H_{int}(t_{1})\nonumber \\{} & {} +(-i)^{2} \int _{t_{i}}^{t_{f}}\!\!\textrm{d}t_{1} \,H_{int}(t_{1}) \int _{t_{i}}^{t_{1}}\textrm{d}t_{2} \, H_{int}(t_{2})+\cdots \nonumber \\ \end{aligned}$$
(13)

where \(H_{int}=\int \!\!\textrm{d}^3 {{{\textbf {x}}}}\,\mathcal {H}_{int}(x)\) is the interaction Hamiltonian.

Notice that we look at the time evolution operator and not at the S-matrix. This is because the phenomenon of flavor oscillations can only be described at finite time. This amounts to say that flavor neutrino states do not exist as asymptotically stable states. As it will be clear from the various examples below, the limits \(t_i \rightarrow - \infty \) and \(t_f \rightarrow + \infty \) forbid the flavor-changing processes under study. At the same time, such a limit guarantees strict energy conservation. This is in agreement with the flavor-energy uncertainty relation derived in Ref. [34] and it is analogous to what happens for unstable particles [27, 28, 35,36,37,38] (see also [39, 40], where the importance of finite-time QFT in the study of decay has been emphasized). As a matter of fact, both the decay of unstable particles [41] and neutrino oscillations [42] can be viewed in terms of the time-energy uncertainty relations.

In the following we will first study the case of 0+1D QFT (that is, QM), and a 3+1D scalar model. This preliminary analysis permits to grasp the main features of the problem, without the complication of dealing with spinors.

3 A quantum mechanics toy model of flavor mixing

Let us consider the quantum mechanical problem of two interacting harmonic oscillators with bare frequencies \(\omega _{A,B}\). We treat this problem as a \(0+1D\) field theory described by the Lagrangian

$$\begin{aligned} L= & {} \frac{1}{2}\left( \frac{dx_{A}}{dt}\right) ^{2}-\frac{\omega _{A}^{2}}{ 2}x_{A}^{2}+\frac{1}{2}\left( \frac{dx_{B}}{dt}\right) ^{2}\nonumber \\{} & {} -\frac{\omega _{B}^{2}}{ 2}y^{2}-\omega _{AB}^{2}x_{A}x_{B} \,. \end{aligned}$$
(14)

In agreement with the previous discussion, we regard the term \(L _{int}=\omega _{AB}^{2}x_{A}x_{B}\) as an interaction, where \(\omega _{AB}^{2}\) (with dimension Energy\(^{2}\)) plays the role of the coupling constant. Hence, the fields in the interaction picture take the form:

$$\begin{aligned} x_{A}(t)= & {} \frac{1}{\sqrt{2\omega _{A}}}\left( a_{A}e^{-i\omega _{A}t}+a_{A}^{\dagger }e^{i\omega _{A}t}\right) \text {,} \end{aligned}$$
(15)
$$\begin{aligned} x_{B}(t)= & {} \frac{1}{\sqrt{2\omega _{B}}}\left( a_{B}e^{-i\omega _{B}t}+a_{B}^{\dagger }e^{i\omega _{B}t}\right) \text {,} \end{aligned}$$
(16)

in which the creation and annihilation operators (with usual commutation relations \([a_{A},a_{A}^{\dagger }]=[a_{B},a_{B}^{\dagger }]=1\) and zero otherwise) have been introduced.

We can safely perform calculations by means of the formula (13), taking \(H_{int}(t)=\omega _{AB}^{2}x_{A}(t)x_{B}(t)\). As initial state \(t_{i}\), we consider an excitation along the A-direction: \(\vert A\rangle =a_{A}^{\dagger }\vert 0\rangle \). We then evaluate the probability that the state has changed at the time \( t_{f}>t_{i},\) a situation that roughly speaking corresponds to a decay of the initial state. The first possible transition is the mixing \(\vert A\rangle =a_{A}^{\dagger }\vert 0\rangle \rightarrow a_{B}^{\dagger }\vert 0\rangle =\vert B\rangle \) driven by the interaction term. The corresponding amplitude reads:

$$\begin{aligned}{} & {} \langle B\vert U(t_{f},t_{i})\vert A\rangle =\langle 0\vert a_{B}U(t_{f},t_{i})a_{A}^{\dagger }\vert 0\rangle \nonumber \\= & {} -i\frac{\omega _{AB}^{2}}{\sqrt{2\omega _{A}}\sqrt{2\omega _{B}}} \int _{t_{i}}^{t_{f}}dt_{1}e^{-i(\omega _{A}-\omega _{B})t_{1}} \nonumber \\= & {} \frac{\omega _{AB}^{2}}{\sqrt{2\omega _{A}}\sqrt{2\omega _{B}}}\frac{ e^{-i(\omega _{A}-\omega _{B})t_{f}}-e^{-i(\omega _{A}-\omega _{B})t_{i}}}{(\omega _{A}-\omega _{B})}\text {.} \end{aligned}$$
(17)

Hence, the probability for this “transition” to happen is:

$$\begin{aligned} \mathcal {P}_{A\rightarrow B}(\Delta t)=\frac{\omega _{AB}^{4}}{\omega _{A}\omega _{B}}\frac{\sin ^{2}\left[ \frac{(\omega _{A}-\omega _{B})\Delta t}{2}\right] }{(\omega _{A}-\omega _{B})^{2}},\end{aligned}$$
(18)

with

$$\begin{aligned} \Delta t=t_{f}-t_{i} . \end{aligned}$$

The formula includes an oscillation whose frequency is proportional to the frequency differences, that we shall call the “low frequency” term. Note, for short times \(\mathcal {P}_{A\rightarrow B}(\Delta t)\simeq \frac{\omega _{AB}^{4}\Delta t^{2}}{ 4\omega _{A}\omega _{B}}.\)

There is, however, at first order another possible transition: \( a_{A}^{\dagger }\vert 0\rangle \rightarrow \frac{\left( a_{A}^{\dagger }\right) ^{2}}{\sqrt{2}}a_{B}^{\dagger }\left| 0\right\rangle ,\) that is a single excitation along A converts into AAB. The corresponding amplitude reads:

$$\begin{aligned}{} & {} \frac{1}{\sqrt{2}}\langle 0\vert a_{B}a_{A}^{2} \, U(t_{f},t_{i})\, a_{A}^{\dagger }\vert 0\rangle \nonumber \\ {}{} & {} \quad =-i \frac{\sqrt{2}\omega _{AB}^{2}}{\sqrt{2\omega _{A}}\sqrt{2\omega _{B}}} \int _{t_{i}}^{t_{f}}dt_{1}e^{-i(\omega _{A}+\omega _{B})t} \nonumber \\{} & {} \quad =\frac{\sqrt{2}\omega _{AB}^{2}}{\sqrt{2\omega _{1}}\sqrt{2\omega _{2}}}\frac{ e^{-i(\omega _{A}+\omega _{B})t_{f}}-e^{-i(\omega _{A}+\omega _{B})t_{i}}}{(\omega _{A}+\omega _{B})}\text {,} \end{aligned}$$
(19)

hence

$$\begin{aligned} \mathcal {P}_{A\rightarrow AAB}(\Delta t)=\frac{2\omega _{AB}^{4}}{\omega _{A}\omega _{B}}\frac{\sin ^{2}\left[ \frac{(\omega _{A}+\omega _{B})\Delta t}{2}\right] }{(\omega _{A}+\omega _{B})^{2}}\text {,} \end{aligned}$$
(20)

which involves the sum of the frequencies and is denoted as the ‘high frequency’ term. For short times, \(P^{A\rightarrow AAB}(\Delta t)\simeq \frac{ \omega _{AB}^{4}t^{2}}{2\omega _{A}\omega _{B}}\).

Summarizing, the total transition probability (in other words, the A transition probability (loosely speaking its decay probability) is given as the sum of both terms.

$$\begin{aligned}{} & {} \mathcal {P}_{D}^{A}(\Delta t) = \mathcal {P}_{A\rightarrow B}(\Delta t)+\mathcal {P}_{A\rightarrow AAB}(\Delta t)\ \nonumber \\{} & {} \quad =\ \frac{\omega _{AB}^{4}}{ \omega _{A}\omega _{B}}\left[ \frac{\sin ^{2}\left[ \frac{(\omega _{A}-\omega _{B})\Delta t}{2} \right] }{(\omega _{A}-\omega _{B})^{2}}+2\frac{\sin ^{2}\left[ \frac{ (\omega _{A}+\omega _{B})\Delta t}{2}\right] }{(\omega _{A}+\omega _{B})^{2}}\right] \,.\nonumber \\ \end{aligned}$$
(21)

For short times, \(\mathcal {P}_{D}^{A}(\Delta t)\simeq \frac{3\omega _{AB}^{4}\Delta t^{2}}{4\omega _{A}\omega _{B}}.\)

Similarly, one can easily calculate within the same framework the survival probability. To this end we need to evaluate

$$\begin{aligned} \langle A\vert U(t_{f},t_{i})\vert A\rangle =\langle 0\vert a_{A}U(t_{f},t_{i})a_{A}^{\dagger }\vert 0\rangle \,. \end{aligned}$$
(22)

Up to the second order we get:

$$\begin{aligned}{} & {} \langle A\vert U(t_{f},t_{i})\vert A\rangle =1\nonumber \\{} & {} \quad -i \, \mathcal {T}\langle 0\vert a_{A}\int _{t_{i}}^{t_{f}}dt_{1}H_{int}(t_{1})\int _{t_{i}}^{t_{1}}dt_{2}H_{int}(t_{2})a_{A}^{\dagger }\vert 0\rangle \,. \end{aligned}$$
(23)

Upon using the equality

$$\begin{aligned}{} & {} \langle 0\vert a_{A}H_{int}(t_{1})H_{int}(t_{2})a_{A}^{\dagger }\vert 0\rangle \nonumber \\{} & {} \quad =\frac{2\omega _{AB}^{4}}{4\omega _{A}\omega _{B}} e^{-i(\omega _{A}+\omega _{B})t_{1}}e^{i(\omega _{A}+\omega _{B})t_{2}}\nonumber \\{} & {} \qquad +\frac{\omega _{AB}^{4}}{4\omega _{A}\omega _{B} }e^{i(\omega _{A}-\omega _{B})t_{1}}e^{-i(\omega _{A}-\omega _{B})t_{2}} \,, \end{aligned}$$
(24)

one gets the survival probability of the state \(\vert A\rangle \) as:

$$\begin{aligned}&\mathcal {P}_{A \rightarrow A}(\Delta t)=\left| 1-\frac{\omega _{AB}^{4}}{4\omega _{A}\omega _{B}}\right. \nonumber \\ {}&\quad \left. \left[ 2\frac{t}{ i(\omega _{A}+\omega _{B})}-2\frac{e^{-i(\omega _{A}+\omega _{B})\Delta t}-1}{(\omega _{A}+\omega _{B})^{2}}\right. \right. \nonumber \\&\quad \left. \left. +\frac{t}{-i(\omega _{A}-\omega _{B})}-\frac{e^{i(\omega _{A}-\omega _{B})\Delta t}-1}{(\omega _{A}-\omega _{B})^{2}}\right] \right| ^{2} \end{aligned}$$
(25)
$$\begin{aligned}&=\left| 1-R-iI\right| ^{2}=(1-R-iI)(1-R+iI)\nonumber \\&=1-R+iI-R+R^{2}-iRI-iI-iIR+I^{2} \end{aligned}$$
(26)
$$\begin{aligned}&=1-2R+\cdots \ , \end{aligned}$$
(27)

where R and I are real. In particular:

$$\begin{aligned} R= & {} \frac{\omega _{AB}^{4}}{4\omega _{A}\omega _{B}}\left( \frac{1-\cos \left[ (\omega _{A}-\omega _{B})\Delta t\right] }{(\omega _{A}-\omega _{B})^{2}}\right. \nonumber \\{} & {} \left. +2\frac{1-\cos \left[ (\omega _{A}+\omega _{B})\Delta t\right] }{(\omega _{A}+\omega _{B})^{2}}\right) \nonumber \\= & {} \frac{\omega _{AB}^{4}}{2\omega _{A}\omega _{B}}\left( \frac{\sin ^{2}\left[ \frac{ (\omega _{A}-\omega _{B})\Delta t}{2}\right] }{(\omega _{A}-\omega _{B})^{2}}+2\frac{\sin ^{2}\left[ \frac{(\omega _{A}+\omega _{B})\Delta t}{2}\right] }{(\omega _{A}+\omega _{B})^{2}}\right) \,, \nonumber \\ \end{aligned}$$
(28)

hence

$$\begin{aligned}{} & {} \mathcal {P}_{S}^{A}(\Delta t) = \mathcal {P}_{A \rightarrow A}(\Delta t) \nonumber \\= & {} \ 1-\frac{\omega _{AB}^{4}}{\omega _{A}\omega _{B}}\left( \frac{\sin ^{2}\left[ \frac{(\omega _{A}-\omega _{B})\Delta t}{2}\right] }{(\omega _{A}-\omega _{B})^{2}}+2\frac{\sin ^{2} \left[ \frac{(\omega _{A}+\omega _{B})\Delta t}{2}\right] }{(\omega _{A}+\omega _{B})^{2}}\right) \,, \nonumber \\ \end{aligned}$$
(29)

which leads to (at order \(g^{2}\)):

$$\begin{aligned} \mathcal {P}_{S}^{A}(\Delta t)+\mathcal {P}_{D}^{A}(\Delta t)\ =\ 1 \,, \end{aligned}$$
(30)

for each t,  as it must. For small t,  \(p_{S}^{A}(t)\simeq 1-\frac{ \omega _{AB}^{4}}{\omega _{A}\omega _{B}}\left( 2\frac{\Delta t^{2}}{4}+\frac{\Delta t^{2}}{4} \right) =1-\frac{3\omega _{AB}^{4}\Delta t^{2}}{4\omega _{A}\omega _{B}}.\)

Of course, the present problem can be also solved by introducing the rotation

$$\begin{aligned} \left( \begin{array}{c} x_{1} \\ x_{2} \end{array} \right) =\left( \begin{array}{cc} \cos \theta &{} -\sin \theta \\ \sin \theta &{} \cos \theta \end{array} \right) \left( \begin{array}{c} x_{A} \\ x_{B} \end{array} \right) \end{aligned}$$
(31)

with

$$\begin{aligned} \theta =\frac{1}{2}\arctan \frac{2\omega _{AB}^{2}}{\omega _{B}^{2}-\omega _{A}^{2}} \,, \end{aligned}$$
(32)

and

$$\begin{aligned} \omega _{1}^{2}&=\omega _{A}^{2}\cos ^{2}\theta +\omega _{B}^{2}\sin ^{2}\theta -\omega _{AB}^{2}\sin (2\theta )\text { }, \end{aligned}$$
(33)
$$\begin{aligned} \omega _{2}^{2}&=\omega _{A}^{2}\sin ^{2}\theta +\omega _{B}^{2}\cos ^{2}\theta +\omega _{AB}^{2}\sin (2\theta )\text { }\,, \end{aligned}$$
(34)
$$\begin{aligned} \omega _{A}^{2}&=\omega _{1}^{2}\cos ^{2}\theta +\omega _{2}^{2}\sin ^{2}\theta \text { }\,, \end{aligned}$$
(35)
$$\begin{aligned} \omega _{B}^{2}&=\omega _{1}^{2}\sin ^{2}\theta +\omega _{2}^{2}\cos ^{2}\theta \text { }\, \text {.} \end{aligned}$$
(36)

The position operators become

$$\begin{aligned} x_{1}(t)= & {} \frac{1}{\sqrt{2\omega _{1}}}\left( a_{1}e^{-i\omega _{1}t}+a_{1}^{\dagger }e^{i\omega _{1}t}\right) \,, \end{aligned}$$
(37)
$$\begin{aligned} x_{2}(t)= & {} \frac{1}{\sqrt{2\omega _{2}}}\left( a_{2}e^{-i\omega _{2}t}+a_{2}^{\dagger }e^{i\omega _{2}t}\right) \,. \end{aligned}$$
(38)

Upon denoting \(\left| \Omega \right\rangle \) as the vacuum of the full Hamiltonian (\(a_{1}\left| \Omega \right\rangle =a_{2}\left| \Omega \right\rangle =0\)), one may also consider the state

$$\begin{aligned} \left| a\right\rangle =\cos \theta a_{1}^{\dagger }\left| \Omega \right\rangle +\sin \theta a_{2}^{\dagger }\left| \Omega \right\rangle \,, \end{aligned}$$
(39)

yet it is clear that \(\left| a\right\rangle \ne \left| A\right\rangle =a_{A}^{\dagger }\left| a\right\rangle \),

In terms of \(\left| a\right\rangle \), the survival probability takes the form:

$$\begin{aligned} \mathcal {P}_{S}^{a}(\Delta t) \ = \ 1-\sin ^{2}2\theta \sin ^{2}\left[ \frac{(\omega _{1}-\omega _{2})\Delta t}{ 2}\right] \, \text {.} \end{aligned}$$
(40)

In the limit of small \(\theta ,\) the previous expression is approximated by:

$$\begin{aligned} \mathcal {P}_{S}^{a}(\Delta t) \ {}\simeq & {} \ 1-\frac{4\omega _{AB}^{4}}{\left( \omega _{B}^{2}-\omega _{A}^{2}\right) ^{2}} \sin ^{2}\left[ \frac{(\omega _{A}-\omega _{B})\Delta t}{2}\right] \nonumber \\= & {} 1-\frac{4\omega _{AB}^{4}}{ \left( \omega _{B}+\omega _{A}\right) ^{2}}\frac{\sin ^{2}\left[ \frac{ (\omega _{A}-\omega _{B})\Delta t}{2}\right] }{\left( \omega _{B}-\omega _{A}\right) ^{2}} \, \text {.} \end{aligned}$$
(41)

We then realize that the functions \(\mathcal {P}_{S}^{a}(\Delta t) \) and \(\mathcal {P}_{S}^{A}(\Delta t) \) are different in various ways. First, the expression \(\mathcal {P}_{S}^{a}(\Delta t) \) contains only the low frequency term but not the high frequency one. Second, the ratio of the coefficients in front of the terms with frequency \(\left( \omega _{A}-\omega _{B}\right) \) reads \(4\omega _{A}\omega _{B}/\left( \omega _{A}+\omega _{B}\right) ^{2}\), which is in general different from unity (it approaches for it in the limit of equal bare masses). This discrepancy is due to the fact that the states \( \left| a\right\rangle \) and \(\left| A\right\rangle \) are different, thus they have different survival probabilities. In the framework of QM, one may “engineer” both initial states. In QFT it is different, and it is not a priori clear to what the field \(x_{A}\) corresponds to.

4 Scalar field mixing in the interaction picture

We now move from QM to QFT. To this end, we investigate the mixing for scalar fields in the interaction picture.

Let us consider two fields \(\phi _{A}=\phi _{A}(t,{\textbf{x}})\) and \(\phi _{B}=\phi _{B}(t, {\textbf{x}})\) that correspond to our flavor bare states A and B,  whose Lagrangian density is given by

$$\begin{aligned} \mathcal {L}= & {} \frac{1}{2}\left( \partial _{\alpha }\phi _{A}\right) ^{2}-\frac{ m_{A}^{2}}{2}\phi _{A}^{2}+\frac{1}{2}\left( \partial _{\alpha }\phi _{B}\right) ^{2}\nonumber \\{} & {} -\frac{m_{B}^{2}}{2}\phi _{B}^{2}-m_{AB}^{2}\phi _{A}\phi _{B}\text {.} \end{aligned}$$
(42)

The Hamiltonian density is

$$\begin{aligned} \mathcal {H}= & {} \frac{\pi _{A}}{2}+\frac{\pi _{B}}{2}+\frac{1}{2}\left( \nabla \phi _{A}\right) ^{2}\nonumber \\{} & {} +\frac{1}{2}\left( \nabla \phi _{B}\right) ^{2}+\frac{m_{A}^{2}}{2}\phi _{A}^{2}+\frac{m_{A}^{2}}{2} \phi _{B}^{2}+m_{AB}^{2}\phi _{A}\phi _{B} \,, \end{aligned}$$
(43)

with \(\pi _{A}=\partial _{t}\phi _{A}\) and \(\pi _{B}=\partial _{t}\phi _{B}.\) We regard the mixing term as a perturbation, thus:

$$\begin{aligned} \mathcal {H}_{0}= & {} \frac{\pi _{A}}{2}+\frac{\pi _{B}}{2}+\frac{1}{2}\left( \nabla \phi _{A}\right) ^{2}+\frac{1}{2}\left( \nabla \phi _{B}\right) ^{2}\nonumber \\{} & {} +\frac{m_{A}^{2}}{2}\phi _{A}^{2}+\frac{m_{A}^{2}}{2}\phi _{B}^{2}, \nonumber \\ \mathcal {H}_{int}= & {} m_{AB}^{2}\phi _{A}\phi _{B} \, \text {.} \end{aligned}$$
(44)

Upon quantizing the system, in the interaction picture we get for the field A:

$$\begin{aligned} \phi _{A}(x)&=\phi _{A}(t,{\textbf{x}})\nonumber \\&=\frac{1}{\sqrt{V}}\sum _{{\textbf {k}}=2\pi {\textbf{n}}/L}\frac{1}{\sqrt{2\omega _{{{\textbf {k}}},A}}}\left( a_{{{\textbf {k}}},A }e^{-ikx}+a_{{{\textbf {k}}},A}^{\dagger }e^{ikx}\right) \, , \end{aligned}$$
(45)
$$\begin{aligned} \pi _{A}(x)&=\pi _{A}(t,{\textbf{x}})\nonumber \\&=\frac{-i}{\sqrt{V}}\sum _{{\textbf {k}} =2\pi {\textbf{n}}/L}\sqrt{\frac{\omega _{{{\textbf {k}}},A}}{2}}\left( a_{{{\textbf {k}}}, A}e^{-ikx}-a_{{{\textbf {k}}},A}^{\dagger }e^{ikx}\right) \, , \end{aligned}$$
(46)

with \(k^{0}=\omega _{{{\textbf {k}}},A}=\sqrt{{\textbf {k}}^{2}+m_{A}^{2}}.\) The commutation relation \(\left[ \phi _{A}(t,{\textbf{x}}),\pi _{A}(t,{\textbf{y}}) \right] =\frac{1}{V}\sum _{{\textbf {k}}}e^{i\mathbf {k\cdot (x-y)}}=i\delta _V( \mathbf {x-y})\) implies that \(\left[ a_{{{\textbf {k}}},A},a_{{{\textbf {p}}},A} ^{\dagger }\right] =\delta _{{{\textbf {k}}},{{\textbf {p}}}}\), zero otherwise. Analogous expressions hold for \(\phi _{B}(x)\) and \(\pi _{B}(x).\)

The interacting Hamiltonian (in the interaction picture) reads:

(47)

Next, we define the “flavor state” A with three-momentum \({\textbf {p}}\) as:

$$\begin{aligned} \left| A,{\textbf {p}}\right\rangle =a_{{{\textbf {p}}},A}^{\dagger }\left| 0\right\rangle \text {.} \end{aligned}$$
(48)

Assuming that such a state is created at \(t=0\), we evaluate the probability that it has transformed into a different state at the time \(t>0\) or, conversely, that it has not changed.

For the case of the transition into a different state, let us first calculate the probability amplitude for the transition \(\left| A,{\textbf{p}}\right\rangle \rightarrow \left| B,{\textbf {k}}\right\rangle \):

(49)

The probability that a particle A with momentum \({\textbf {p}}\) converts into a particle B is obtained upon summing over the density of final states \( \sum _{{\textbf {k}}}\):

(50)

The result is finite and well behaved. Note, for short time \(\mathcal {P}_{A\rightarrow B}({\textbf {p}}; \Delta t)\simeq \frac{m_{AB}^{4}}{4\omega _{{{\textbf {p}}},A} \omega _{{{\textbf {p}}},B}}t^{2}.\)

Yet, other transitions are possible. Namely, we may have the transition \( A\rightarrow AAB\) of the type

$$\begin{aligned} \left| A,{\textbf {p}}\right\rangle \rightarrow \left| A,{\textbf {k}} _{1}\right\rangle \left| A,{\textbf {k}}_{2}\right\rangle \left| B, {\textbf {k}}_{3}\right\rangle =a_{A,{\textbf {k}}_{1}}^{\dagger }a_{A,{\textbf {k}} _{2}}^{\dag }a_{B,{\textbf {k}}_{3}}^{\dag }\left| 0\right\rangle \text {,} \end{aligned}$$
(51)

where the two emitted A particles have different momentum, \({\textbf {k}} _{1}\ne {\textbf {k}}_{2}.\) The corresponding amplitude of this process reads

$$\begin{aligned}{} & {} \mathcal {A}_{A\rightarrow AAB}^{{\textbf {k}}_{1}\ne {\textbf {k}}_{2}}({\textbf {p}},{\textbf {k}}_{1},{\textbf {k}}_{2},{\textbf {k}}_{3};t_i,t_f)\nonumber \\{} & {} \quad =-i\int _{t_{i}}^{t_{f}}dt_{1}\langle 0\vert a_{A,{\textbf {k}} _{1}}a_{A,{\textbf {k}}_{2}}a_{B,{\textbf {k}}_{3}}H_{int}(t_{1})a_{A,{\textbf {p}} }^{\dagger }\vert 0\rangle \text {.} \end{aligned}$$
(52)

After an explicit calculation up to first order, its squared modulus turns out to be:

$$\begin{aligned}{} & {} \vert \mathcal {A}_{A\rightarrow AAB}^{{\textbf {k}}_{1}\ne {\textbf {k}} _{2}}({\textbf {p}},{\textbf {k}}_{1},{\textbf {k}}_{2},{\textbf {k}} _{3},t_i,t_f)\vert ^{2}\nonumber \\{} & {} \quad =\frac{m_{AB}^{4}}{\omega _{{{\textbf {k}}}_3,A} \omega _{{{\textbf {k}}}_3,B}}\frac{\sin ^{2}\Big [ \frac{\left( \omega _{{\textbf {k}}_{3},A} +\omega _{{{\textbf {k}}}_3,B}\right) \Delta t}{2}\Big ] }{ \left( \omega _{{\textbf {k}}_{3},A} +\omega _{{{\textbf {k}}}_3,B}\right) ^{2}}\times \nonumber \\ {}{} & {} \quad \times \left( \delta _{{\textbf {k}}_{1},{\textbf {p}}}\delta _{{\textbf {k}}_{2},-{\textbf {k}} _{3}}+\delta _{{\textbf {k}}_{1},-{\textbf {k}}_{3}}\delta _{{\textbf {k}}_{2}, {\textbf {p}}}\right) \text {.} \end{aligned}$$
(53)

Next, one needs to sum over final states \({\textbf {k}}_{1},{\textbf {k}}_{2}, {\textbf {k}}_{3}.\) leading to the probability:

$$\begin{aligned}{} & {} \mathcal {P}^{{{\textbf {k}}}_1 \ne {{\textbf {k}}}_2}_{A\rightarrow AAB}({\textbf {p}};\Delta t) \nonumber \\ {}{} & {} \quad =\text { }\frac{1}{2}\sum _{ {\textbf {k}}_{1},{\textbf {k}}_{2},{\textbf {k}}_{3}}\left| \mathcal {A}_{A\rightarrow AAB}^{{\textbf {k}}_{1}\ne {\textbf {k}}_{2}}({\textbf {p}},{\textbf {k}}_{1},{\textbf {k}}_{2},{\textbf {k}}_{3},t_i,t_f)\right| ^{2} \nonumber \\= & {} \sum _{{\textbf {k}}_{3}}\frac{m_{AB}^{4}}{\omega _{{{\textbf {k}}}_3,A}\omega _{{{\textbf {k}}}_3,B}}\frac{\sin ^{2}\Big [ \frac{\left( \omega _{{{\textbf {k}}}_3,A} +\omega _{{{\textbf {k}}}_3,B}\right) \Delta t}{2}\Big ] }{\left( \omega _{{\textbf {k}}_{3},A} +\omega _{{{\textbf {k}}}_3,B}\right) ^{2}}\nonumber \\{} & {} -\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}} \,, \end{aligned}$$
(54)

where the factor 1/2 in front of the sum takes into account that the two A in the final state are identical bosons. The subtracted term in the last equation is due to the condition \({\textbf {k}}_{2}\ne {\textbf {k}}_{1}.\) Note, the sum term diverges, thus a certain cutoff is implicitly introduced so to keep the intermediate results finite (which is then sent to infinity at the very end of the calculation). When the volume is sufficiently large, the previous expression becomes

$$\begin{aligned}{} & {} \mathcal {P}^{{{\textbf {k}}}_1 \ne {{\textbf {k}}}_2}_{A\rightarrow AAB}({{\textbf {p}}};\Delta t)\nonumber \\{} & {} \quad =V\int \!\! \frac{\textrm{d}^{3}{{\textbf {k}}}_{3}}{(2\pi )^{3}} \, \frac{m_{AB}^{4}}{\omega _{{{\textbf {k}}}_3,A}\omega _{{{\textbf {k}}}_3,B}}\frac{\sin ^{2}\Big [ \frac{\left( \omega _{{{\textbf {k}}}_3,A} +\omega _{{{\textbf {k}}}_3,B}\right) \Delta t}{2}\Big ] }{\left( \omega _{{\textbf {k}}_{3},A} +\omega _{{{\textbf {k}}}_3,B}\right) ^{2}}\nonumber \\{} & {} \qquad -\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}} \,, \end{aligned}$$
(55)

where, again, a cutoff is implicit in the integral over \(k_{3}\). The term proportional to V is a typical vacuum term that needs to be subtracted. However, the second term in Eq. (55) needs to be kept, see below.

The last possible transition is the case in which \({\textbf {k}}_{2}= {\textbf {k}}_{1},\) thus

$$\begin{aligned} \vert A,{\textbf {p}}\rangle \rightarrow \frac{1}{\sqrt{2}}a_{{{\textbf {k}}}_1, A}^{\dagger }a_{{\textbf {k}}_{1},A}^{\dag }a_{{\textbf {k}}_{3},B }^{\dag }\vert 0\rangle \text {.} \end{aligned}$$
(56)

The amplitude at first order is

$$\begin{aligned}{} & {} \mathcal {A}_{A\rightarrow AAB}^{{\textbf {k}}_{1}={\textbf {k}}_{2}}({\textbf {p}}, {\textbf {k}}_{1},{\textbf {k}}_{3};t_i,t_f)\nonumber \\{} & {} \quad =-\frac{i}{\sqrt{2}} \int _{t_{i}}^{t_{f}}dt_{1}\langle 0\vert a_{{\textbf {k}}_{1},A}a_{ {\textbf {k}}_{1},A}a_{{\textbf {k}}_{3},B}H_{int}(t_{1})a_{{\textbf {p}},A}^{\dagger }\vert 0\rangle , \end{aligned}$$
(57)

whose squared modulus is

$$\begin{aligned}{} & {} \vert \mathcal {A}_{A\rightarrow AAB}^{{\textbf {k}}_{1}={\textbf {k}}_{2}}( {\textbf {p}},{\textbf {k}}_{1},{\textbf {k}}_{3};t_i,t_f)\vert ^{2}\nonumber \\{} & {} \quad =\delta _{ {\textbf {k}}_{1},{\textbf {p}}}\delta _{{\textbf {k}}_{1},\mathbf {-k}_{3}}2\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}} \,. \end{aligned}$$
(58)

Upon summing over the final momenta \({\textbf {k}}_{1},{\textbf {k}}_{3}\):

$$\begin{aligned} \mathcal {P}^{{{\textbf {k}}}_1 = {{\textbf {k}}}_2}_{A\rightarrow AAB}({\textbf {p}},\Delta t)= & {} \sum _{{\textbf {k}}_{1}, {\textbf {k}}_{3}}\left| \mathcal {A}_{A\rightarrow AAB}^{{\textbf {k}}_{1}= {\textbf {k}}_{2}}({\textbf {p}},{\textbf {k}}_{1},{\textbf {k}}_{3};t_i,t_f)\right| ^{2}\nonumber \\= & {} 2\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}}\text {.} \end{aligned}$$
(59)

Note, the factor 2 appears just as in the QM toy model of Sect. 3.

Putting all the pieces together, the total probability of A going into something else (thus, a decay probability) is:

$$\begin{aligned}{} & {} \mathcal {P}_{D}^{A}({\textbf {p}}; \Delta t) =\mathcal {P}_{A\rightarrow B}({\textbf {p}}; \Delta t) +\mathcal {P}^{{{\textbf {k}}}_1 \ne {{\textbf {k}}}_2}_{A\rightarrow AAB}({\textbf {p}}; \Delta t) \nonumber \\{} & {} \quad +\mathcal {P}^{{{\textbf {k}}}_1 = {{\textbf {k}}}_2}_{A\rightarrow AAB}({\textbf {p}}; \Delta t) \nonumber \\{} & {} =\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B} \right) ^{2}}\nonumber \\{} & {} \quad +\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}}\nonumber \\{} & {} \quad +V\int \!\! \frac{\textrm{d}^{3}{{\textbf {k}}}_{3}}{(2\pi )^{3}}\frac{m_{AB}^{4}}{\omega _{{{\textbf {k}}}_3,A} \omega _{{{\textbf {k}}}_3,B}}\frac{\sin ^{2}\Big [ \frac{\left( \omega _{{{\textbf {k}}}_3,A}+\omega _{{{\textbf {k}}}_3,B}\right) \Delta t}{2}\Big ] }{ \left( \omega _{{{\textbf {k}}}_3, A}+\omega _{{{\textbf {k}}}_3,B}\right) ^{2}} \text {.} \end{aligned}$$
(60)

The factor 2 of \(\mathcal {P}^{{{\textbf {k}}}_1 = {{\textbf {k}}}_2}_{A\rightarrow AAB}({\textbf {p}}; \Delta t)\) combines with the factor \(-1\) in \(\mathcal {P}^{{{\textbf {k}}}_1 \ne {{\textbf {k}}}_2}_{A\rightarrow AAB}({\textbf {p}}; \Delta t)\) in order to give the same factor in front of the high-frequency term. Finally, the term proportional to V, being a vacuum term, is subtracted. Then the probability that the oscillation takes place up to second order is:

$$\begin{aligned} \mathcal {P}_{D}^{A}({\textbf {p}};\Delta t)= & {} \frac{m_{AB}^{4}}{\omega _{{{\textbf {p}}},A} \omega _{{{\textbf {p}}},B}}\left( \frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B} \right) ^{2}}\right. \nonumber \\{} & {} \left. +\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}}\right) \text {.} \end{aligned}$$
(61)

It is important to verify the correctness of the previous result. Just as in the QM toy model, one needs to check the survival probability of the state \( \left| A,{\textbf {p}}\right\rangle \). To this end, we calculate the probability of the transition

$$\begin{aligned} \left| A,{\textbf {p}}\right\rangle \rightarrow \left| A,{\textbf {k}} \right\rangle \,, \end{aligned}$$
(62)

and then sum over \(\mathbf {k.}\) The corresponding amplitude, up to second order, reads:

$$\begin{aligned}{} & {} \mathcal {A}_{A\rightarrow A}({\textbf {p}},{\textbf {k}};t_i,t_f)=\langle 0\vert a_{{\textbf {k}},A}^{\dagger }U(t_{f},t_{i})a_{{\textbf {p}},A}^{\dagger }\vert 0\rangle \nonumber \\{} & {} =\delta _{\textbf{k,p}}+(-i)^{2}\int _{t_{i}}^{t_{f}}dt_{1}\int _{t_{i}}^{t_{1}}dt_{2}\langle 0\vert a_{{\textbf {k}},A} H_{int}(t_{1})H_{int}(t_{2})a_{{\textbf {p}},A }^{\dagger }\vert 0\rangle \text {.}\nonumber \\ \end{aligned}$$
(63)

Its modulus square takes the form:

$$\begin{aligned}{} & {} \left| \mathcal {A}_{A\rightarrow A}({\textbf {p}},{\textbf {k}};t_i,t_f)\right| ^{2} \nonumber \\ {}{} & {} =\delta _{\textbf{pk}}\left( 1-\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B} \right) ^{2}}\right. \nonumber \\{} & {} \quad \left. -\frac{m_{AB}^{4} }{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2}\left[ \frac{ \left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2} \right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) ^{2}} \right) \nonumber \\{} & {} \quad -\delta _{\textbf{pk}}\sum _{{\textbf {q}}_{1}}\frac{m_{AB}^{4}}{\omega _{{\textbf {q}}_{1},A} 2\omega _{{\textbf {q}}_{1},B}}\frac{\sin ^{2}\Big [ \frac{ \left( \omega _{{\textbf {q}}_{1},A} +\omega _{{\textbf {q}}_{1},B}\right) \Delta t }{2}\Big ] }{\left( \omega _{{\textbf {q}}_{1},A} +\omega _{{\textbf {q}}_{1},B} \right) ^{2}}+\cdots \, \end{aligned}$$
(64)

where dots refer to higher order terms.

Then, upon summing over the final three-momentum \({\textbf {k}}\) and taking the large volume limit, the survival probability reads (up to second order):

$$\begin{aligned}{} & {} \mathcal {P}^A_{S}({\textbf {p}};\Delta t) = \mathcal {P}_{A \rightarrow A}({\textbf {p;}}\Delta t)\nonumber \\ {}{} & {} =\sum _{{\textbf {k}}}\left| \mathcal {A} _{A\rightarrow A}({\textbf {p}},{\textbf {k}};t_i,t_f)\right| ^{2} \nonumber \\{} & {} =1-\frac{ m_{AB}^{4}}{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B} \right) ^{2}}\nonumber \\{} & {} -\frac{m_{AB}^{4} }{\omega _{{{\textbf {p}}},A}\omega _{{{\textbf {p}}},B}}\frac{\sin ^{2}\left[ \frac{ \left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2} \right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) ^{2}} \nonumber \\{} & {} -V\int \!\! \frac{\textrm{d}^{3}{{\textbf {q}}}_1}{(2\pi )^{3}} \, \frac{m_{AB}^{4}}{\omega _{{\textbf {q}}_{1},A} 2\omega _{{\textbf {q}}_{1},B}}\frac{\sin ^{2}\Big [ \frac{ \left( \omega _{{\textbf {q}}_{1},A} +\omega _{{\textbf {q}}_{1},B}\right) \Delta t }{2}\Big ] }{\left( \omega _{{\textbf {q}}_{1},A} +\omega _{{\textbf {q}}_{1},B} \right) ^{2}} \,, \end{aligned}$$
(65)

where the latter term corresponds, in diagrammatic term, to the disconnected vacuum diagram with an AB loop. This term coincides exactly with the one obtained previously. Upon subtracting this term, we find:

$$\begin{aligned} \mathcal {P}^A_{S}({\textbf {p}};\Delta t)= & {} \ 1- \frac{m_{AB}^{4}}{\omega _{{{\textbf {p}}},A} \omega _{{{\textbf {p}}},B}}\left( \frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} -\omega _{{{\textbf {p}}},B} \right) ^{2}}\right. \nonumber \\{} & {} \left. +\frac{\sin ^{2} \left[ \frac{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B}\right) \Delta t}{2}\right] }{\left( \omega _{{{\textbf {p}}},A} +\omega _{{{\textbf {p}}},B} \right) ^{2}}\right) \,, \end{aligned}$$
(66)

with

$$\begin{aligned} \mathcal {P}^A_{D}({\textbf {p}};\Delta t) +\mathcal {P}^A_{S}({\textbf {p}};\Delta t) \ =\ 1 \,, \end{aligned}$$
(67)

as it must.

We thus obtain the probability that the flavor A oscillates (or does not oscillate) as the sum of two distinct term involving the low-frequency and the high-frequency term, where the frequencies \(\omega _{{{\textbf {p}}},A}\) and \(\omega _{{{\textbf {p}}},B}\) depend on the chosen momentum \(\mathbf {p.}\) Note, the structure of the solution is very similar to the QM toy model besides the factor 2 of the high-frequency term. It turns out that only 1 survives the process of renormalization.

Also in the scalar QFT case one may introduce a suitable diagonalization of the fields

$$\begin{aligned} \left( \begin{array}{c} \phi _{1} \\ \phi _{2} \end{array} \right) =\left( \begin{array}{cc} \cos \theta &{} -\sin \theta \\ \sin \theta &{} \cos \theta \end{array} \right) \left( \begin{array}{c} \phi _{A} \\ \phi _{B} \end{array} \right) \,, \end{aligned}$$
(68)

with

$$\begin{aligned} \theta =\frac{1}{2}\arctan \frac{2m_{AB}^{2}}{m_{B}^{2}-m_{A}^{2}} \,. \end{aligned}$$
(69)

The fields \(\phi _{1}\) and \(\phi _{2}\) contain the annihilation (creation) operators \(a_{1,\textbf{p,}}\) \(a_{2,\textbf{p,}}\) (\(a_{1,{\textbf {p}}}^{\dagger }\),\(a_{2, {\textbf {p}}}^{\dagger })\). If we introduce the “Pontecorvo state”

$$\begin{aligned} \left| a,{\textbf {p}}\right\rangle =\cos \theta a_{1}^{\dagger }\left| \Omega \right\rangle +\sin \theta a_{2}^{\dagger }\left| \Omega \right\rangle \,, \end{aligned}$$
(70)

the corresponding survival probability takes the form:

$$\begin{aligned} \mathcal {P}_S^a({{\textbf {p}}};\Delta t)\ = \ 1-\sin ^{2}2\theta \sin ^{2}\left[ \frac{(\omega _{{{\textbf {p}}},1}-\omega _{{{\textbf {p}}},2})\Delta t}{2}\right] \,, \end{aligned}$$
(71)

with \(\omega _{{{\textbf {p}}},j}=\sqrt{|{\textbf {p}}|^{2}+m_{j}^{2}}\), \(j=1,2\). In the limit of small \(\theta ,\) the previous expression is approximated by:

$$\begin{aligned}{} & {} \mathcal {P}_S^a({{\textbf {p}}};\Delta t) \simeq 1-\frac{4m_{AB}^{4}}{\left( m _{B}^{2}-m_{{{\textbf {A}}}}^{2}\right) ^{2}}\sin ^{2}\left[ \frac{(\omega _{{{\textbf {p}}},A}-\omega _{{{\textbf {p}}},B})\Delta t}{2}\right] \,, \nonumber \\ \end{aligned}$$
(72)

which, just as in the QM toy model, differs from Eq. (66) since the high-frequency term is missing and because the factor in front of the low-frequency one is not the same (but the two expressions degenerate for large \(\left| {\textbf {p}}\right| \)). The difference is expected because \(\left| a,{\textbf {p}}\right\rangle \ne \left| A,{\textbf {p}} \right\rangle .\)

5 Neutrino oscillations in the interaction picture

Let us now deal with the fermion case, which can be naturally applied to neutrino oscillations and it will be then referred as neutrino case.

In the interaction picture \(\nu _\sigma \) (\(\sigma =e,\mu \)), defined by the Lagrangian Eq. (1), can be expanded as free fields, evolving under the action of \(\mathcal {L}_0\):

$$\begin{aligned}{} & {} \nu _{\sigma }(x) = \frac{1}{\sqrt{V}} \sum _{{{\textbf {k}}},r}\, \left[ u_{\textbf{k},\sigma }^{r}(t) \, \alpha _{\textbf{k},\sigma }^{r} + v_{-\textbf{k},\sigma }^{r}(t) \, \beta _{-\textbf{k},\sigma }^{r\dag } \right] e^{i\textbf{k}\cdot \textbf{x}}\,, \nonumber \\ \end{aligned}$$
(73)

with \(u^r_{\textbf{k},\sigma }(t) \,= \, e^{- i \omega _{{{\textbf {k}}},\sigma } t}\, u^r_{\textbf{k},\sigma }\;\), \(\;v^r_{\textbf{k},\sigma }(t) \,= \, e^{ i \omega _{{{\textbf {k}}},\sigma } t}\, v^r_{\textbf{k},\sigma }\), \(\omega _{{{\textbf {k}}},\sigma }=\sqrt{|{{\textbf {k}}}|^2 + m_\sigma ^2}\). Annihilation operators satisfy

$$\begin{aligned} \alpha ^r_{{{\textbf {k}}}, \sigma }|0 \rangle = 0 = \beta _{\textbf{k},\sigma }^{r} |0 \rangle . \end{aligned}$$
(74)

The anticommutation relations are

$$\begin{aligned} \{\alpha ^r_{\textbf{k},\rho }, \alpha ^{s\dag }_{\textbf{q},\sigma }\} = \delta _{{{\textbf {k}}} {{\textbf {q}}}}\delta _{rs}\delta _{\rho \sigma } \quad \,, \quad \{\beta ^r_{\textbf{k},\rho }, \beta ^{s\dag }_{\textbf{q},\sigma }\} = \delta _{{{\textbf {k}}} {{\textbf {q}}}} \delta _{rs}\delta _{\rho \sigma },\nonumber \\ \end{aligned}$$
(75)

and the spinors are normalized so that

$$\begin{aligned} u^{r\dag }_{\textbf{k},\rho } u^{s}_{\textbf{k},\rho } = v^{r\dag }_{\textbf{k},\rho } v^{s}_{\textbf{k},\rho } \ = \ \delta _{rs} \quad , \quad u^{r\dag }_{\textbf{k},\rho } v^{s}_{-\textbf{k},\rho } = 0 . \end{aligned}$$
(76)

As in the previous examples, the idea is to perform the perturbative calculation up the first in \(m_{e \mu }\). The interacting Hamiltonian reads:

$$\begin{aligned} H_{int}(t)= & {} m_{e \mu } \sum _{s,s'=1,2}\sum _{{{\textbf {p}}}} \Big [\beta ^s_{{{\textbf {p}}},\mu }\beta ^{s\dag }_{{{\textbf {p}}},e} \delta _{s s'} W^*_{{\textbf {p}}}(t) \nonumber \\ {}{} & {} +\alpha ^{r\dag }_{{{\textbf {p}}},\mu } \alpha ^r_{{{\textbf {p}}},e} \delta _{s s'} W_{{\textbf {p}}}(t) \nonumber \\ {}{} & {} + \beta ^s_{-{{\textbf {p}}},\mu }\alpha ^{s'}_{e,{{\textbf {p}}}} \left( Y^{s s'}_{{\textbf {p}}}(t)\right) ^*\nonumber \\ {}{} & {} +\alpha ^{s\dag }_{{{\textbf {p}}},\mu }\beta ^{s'\dag }_{-{{\textbf {p}}},e} Y^{s s'}_{{\textbf {p}}}(t)\, + \, e \leftrightarrow \mu \Big ] \,, \end{aligned}$$
(77)

where we defined

$$\begin{aligned} W_{{\textbf {p}}}(t)= & {} {\overline{u}}^s_{{{\textbf {p}}},\mu } u^s_{{{\textbf {p}}},e} e^{i\left( \omega _{{{\textbf {k}}},\mu } -\omega _{{{\textbf {k}}},e}\right) t} \ = \ W_{{\textbf {p}}} \, e^{i\left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) t} \end{aligned}$$
(78)
$$\begin{aligned} Y^{s s'}_{{\textbf {p}}}(t)= & {} \, {\overline{u}}^{s}_{{{\textbf {p}}},\mu } v^{s'}_{-{{\textbf {p}}},e} e^{i\left( \omega _{{{\textbf {k}}},\mu }+\omega _{{{\textbf {k}}},e}\right) t} \ = \ Y^{s s'}_{{\textbf {p}}} e^{i\left( \omega _{{{\textbf {p}}},\mu }+\omega _{{{\textbf {p}}},e}\right) t} \end{aligned}$$
(79)

Explicitly

$$\begin{aligned} W_{{\textbf {p}}}= & {} \sqrt{\frac{\left( \omega _{{{\textbf {p}}},e}+m_{e}\right) \left( \omega _{{{\textbf {p}}},\mu }+m_{\mu }\right) }{4\omega _{{{\textbf {p}}},e} \omega _{{{\textbf {p}}},\mu }}} \nonumber \\{} & {} \times \left( 1-\frac{|{{\textbf {p}}}|^{2}}{(\omega _{{{\textbf {p}}},e}+m_{e})(\omega _{{{\textbf {p}}},\mu }+m_{\mu })}\right) \,, \end{aligned}$$
(80)
$$\begin{aligned}{} & {} Y^{22}_{{{\textbf {p}}}} = -Y^{11}_{{{\textbf {p}}}} \ \nonumber \\{} & {} = \ \frac{p_3}{\sqrt{4 \omega _{{{\textbf {p}}},e}\omega _{{{\textbf {p}}},\mu }}} \left( \sqrt{\frac{\omega _{{{\textbf {p}}},\mu }+m_\mu }{\omega _{{{\textbf {p}}},e}+m_e}}+\sqrt{\frac{\omega _{{{\textbf {p}}},e}+m_e}{\omega _{{{\textbf {p}}},\mu }+m_\mu }}\right) \,, \nonumber \\ \end{aligned}$$
(81)
$$\begin{aligned}{} & {} Y^{12}_{{{\textbf {p}}}} = \left( Y^{21}_{{{\textbf {p}}}} \right) ^* \ \nonumber \\{} & {} = \ -\frac{p_1-i p_2}{\sqrt{4 \omega _{{{\textbf {p}}},e}\omega _{{{\textbf {p}}},\mu }}} \left( \sqrt{\frac{\omega _{{{\textbf {p}}},\mu }+m_\mu }{\omega _{{{\textbf {p}}},e}+m_e}}+\sqrt{\frac{\omega _{{{\textbf {p}}},e}+m_e}{\omega _{{{\textbf {p}}},\mu }+m_\mu }}\right) \,.\nonumber \\ \end{aligned}$$
(82)

The first non-trivial flavor transition process we consider is

$$\begin{aligned} |\nu ^r_{{{\textbf {p}}},e}\rangle \ \rightarrow \ |\nu ^s_{{{\textbf {k}}},\mu }\rangle \,, \qquad |\nu ^r_{{{\textbf {p}}},\sigma }\rangle \equiv \alpha ^{r\dag }_{{{\textbf {p}}},\sigma }|0\rangle \,, \end{aligned}$$
(83)

whose amplitude reads

$$\begin{aligned}{} & {} \mathcal {A}^{rs}_{e \rightarrow \mu }({{\textbf {p}}},{{\textbf {k}}},;t_i,t_f) \ \approx \ - i m_{e \mu }\delta _{r s} \delta _{{{\textbf {k}}},{{\textbf {p}}}} W_{{{\textbf {p}}}}\, \int ^{t_f}_{t_i} \!\! \textrm{d}t \, e^{i\left( \omega _{{{\textbf {k}}},\mu }-\omega _{{{\textbf {p}}},e}\right) t} \nonumber \\{} & {} \quad = m_{e \mu } \, \delta _{r s} \delta _{{{\textbf {k}}},{{\textbf {p}}}} \, \, \left( e^{i\left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) t_f}-e^{i \left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) t_i}\right) \frac{W_{{{\textbf {p}}}}}{\omega _{{{\textbf {k}}},e}-\omega _{{{\textbf {k}}},\mu }} \nonumber \\ {}{} & {} \quad = \ \delta _{r s} \delta _{{{\textbf {k}}},{{\textbf {p}}}} \, \tilde{\mathcal {A}}_{e \rightarrow \mu }({{\textbf {k}}}; t_i,t_f), \end{aligned}$$
(84)

where

$$\begin{aligned} \tilde{\mathcal {A}}_{e \rightarrow \mu }({{\textbf {p}}}; t_i,t_f)= & {} \frac{m_{e \mu } \, W_{{{\textbf {p}}}}}{\omega _{{{\textbf {p}}},e}-\omega _{{{\textbf {p}}},\mu }} \, \left( e^{i\left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) t_f}\right. \nonumber \\{} & {} \quad \left. -e^{i \left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) t_i}\right) \,. \end{aligned}$$
(85)

Similarly as in the boson case (see Eq. (50)), the oscillation probability is computed by summing over the final density of states, now involving the sum over the helicities:

$$\begin{aligned}{} & {} \mathcal {P}_{e \rightarrow \mu }({{\textbf {p}}};\Delta t) = \sum _{{{\textbf {k}}},s} |\mathcal {A}^{rs}_{e \rightarrow \mu }({{\textbf {p}}}, {{\textbf {k}}};t_i,t_f)|^2 \nonumber \\ {}{} & {} = |\tilde{\mathcal {A}}_{e \rightarrow \mu }({{\textbf {p}}}, t_i,t_f)|^2 \nonumber \\{} & {} = W_{{{\textbf {p}}}}^2 \, \frac{2 m^2_{e \mu } }{\left( \omega _{{{\textbf {p}}},e}-\omega _{{{\textbf {p}}},\mu }\right) ^2} \left[ 1-\cos \left[ \left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) \Delta t\right] \right] . \end{aligned}$$
(86)

Another non-trivial process is the decay

$$\begin{aligned} |\nu ^r_{{{\textbf {p}}},e}\rangle \ \rightarrow \ |\nu ^{s_1}_{{{\textbf {k}}}_1,e}\rangle |\nu ^{s_2}_{{{\textbf {k}}}_2,\mu }\rangle |\overline{\nu }^{s_3}_{{{\textbf {k}}}_3,e}\rangle .\nonumber \\ \end{aligned}$$
(87)

The amplitude explicitly reads

$$\begin{aligned}{} & {} \mathcal {A}^{r s_1 s_2 s_3}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {p}}},{{\textbf {k}}}_1,{{\textbf {k}}}_2, {{\textbf {k}}}_3;t_i,t_f) \nonumber \\{} & {} \quad \approx -i \, m_{e \mu } \, \, Y_{{{\textbf {k}}}_2}^{s_3 s_2} \, \delta _{{{\textbf {k}}}_1, {{\textbf {p}}}} \delta _{{{\textbf {k}}}_2, - {{\textbf {k}}}_3}\, \delta _{r s_1} \int ^{t_f}_{t_i} \!\! \textrm{d}t \, e^{-i\left( \omega _{{{\textbf {k}}}_2,\mu }+\omega _{{{\textbf {k}}}_2,e}\right) t} \nonumber \\{} & {} \quad = \ -m_{e \mu } \, \delta _{r s_1} \, \delta _{{{\textbf {k}}}_1, {{\textbf {p}}}} \delta _{{{\textbf {k}}}_2, - {{\textbf {k}}}_3} \, \left( e^{-i\left( \omega _{{{\textbf {k}}}_2,\mu }+\omega _{{{\textbf {k}}}_2,e}\right) t_f}\right. \nonumber \\{} & {} \qquad \left. -e^{-i \left( \omega _{{{\textbf {k}}}_2,\mu }+\omega _{{{\textbf {k}}}_2,e}\right) t_i}\right) \frac{Y_{{{\textbf {k}}}_2}^{s_2 s_3}}{\omega _{{{\textbf {k}}}_2,e}+\omega _{{{\textbf {k}}}_3,\mu }} \nonumber \\{} & {} \quad = \ \delta _{{{\textbf {k}}}_1, {{\textbf {p}}}} \delta _{{{\textbf {k}}}_2, - {{\textbf {k}}}_3} \, \delta _{r s_1} \, \tilde{\mathcal {A}}^{s_2 s_3}_{e \rightarrow e\overline{\mu }\mu }({{\textbf {k}}}_2;t_i,t_f) \,, \end{aligned}$$
(88)

where

$$\begin{aligned}{} & {} \tilde{\mathcal {A}}^{s_2 s_3}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {k}}};t_i,t_f) \nonumber \\{} & {} \quad = \ -\frac{m_{e \mu } \, Y^{s_2 s_3}_{{{\textbf {k}}}}}{\omega _{{{\textbf {k}}},e}+\omega _{{{\textbf {k}}},\mu }} \, \left( e^{-i\left( \omega _{{{\textbf {k}}},\mu }+\omega _{{{\textbf {k}}},e}\right) t_f}-e^{-i \left( \omega _{{{\textbf {k}}},\mu }+\omega _{{{\textbf {k}}},e}\right) t_i}\right) \,.\nonumber \\ \end{aligned}$$
(89)

As done above, we thus find the probability as

$$\begin{aligned}{} & {} \mathcal {P}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {p}}};\Delta t) \nonumber \\{} & {} \quad = \ \sum _{{{\textbf {k}}}_1,{{\textbf {k}}}_2,{{\textbf {k}}}_3} \sum _{s_1,s_2,s_3} |\mathcal {A}^{r s_1 s_2 s_3}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {p}}},{{\textbf {k}}}_1,{{\textbf {k}}}_2, {{\textbf {k}}}_3;t_i,t_f) |^2 \nonumber \\{} & {} \quad = \ \sum _{{{\textbf {k}}}}\sum _{s_2,s_3}|\tilde{\mathcal {A}}^{s_2 s_3}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {k}}};t_i,t_f)|^2\,. \end{aligned}$$
(90)

In the large-V limit

$$\begin{aligned}{} & {} \mathcal {P}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {p}}};\Delta t) \nonumber \\{} & {} \quad = \ V \sum _{s_2,s_3} \, \int \!\! \frac{\textrm{d}^3 {{\textbf {k}}}}{(2 \pi )^3} \, \frac{\left( Y^{s_2 s_3}_{{\textbf {k}}}\right) ^2 }{\left( \omega _{{{\textbf {k}}},e}+\omega _{{{\textbf {k}}},\mu }\right) ^2}\times \nonumber \\ {}{} & {} \qquad \times \sin ^2\left( \frac{\left( \omega _{{{\textbf {k}}},\mu }+\omega _{{{\textbf {k}}},e}\right) \Delta t}{2}\right) \,. \end{aligned}$$
(91)

This is a divergent contribution (vacuum diagram) and it must be subtracted.

Finally, we have the process

$$\begin{aligned} |\nu ^r_{{{\textbf {p}}},e}\rangle \ \rightarrow \ |\nu ^{s_1}_{{{\textbf {k}}}_1,e}\rangle |\nu ^{s_2}_{{{\textbf {k}}}_2,e}\rangle |\overline{\nu }^{s_3}_{{{\textbf {k}}}_3,\mu }\rangle \,, \qquad {{\textbf {k}}}_1 \ne {{\textbf {k}}}_2 \ \vee \ s_1 \ne s_2 \,.\nonumber \\ \end{aligned}$$
(92)

The amplitude explicitly reads

(93)

where \(\tilde{\mathcal {A}}^{s_2 s_3}_{e \rightarrow e e \overline{\mu } }({{\textbf {k}}};t_i,t_f)=\tilde{\mathcal {A}}^{s_2 s_3}_{e \rightarrow e{\overline{e}} \mu }({{\textbf {k}}};t_i,t_f)\). Note this correctly goes to zero when \({{\textbf {k}}}_1={{\textbf {k}}}_2\) and \(s_1=s_2\). We thus find the probability as

(94)

It is clear we can not simply subtract the first piece involving the sum over momenta: this procedure would give a negative probability. This subtle issue can be solved by remembering that asymptotic states have to be representations of momentum. Thus, we must isolate the contribution with \({{\textbf {k}}}= {{\textbf {p}}}\)

$$\begin{aligned}{} & {} \mathcal {P}_{e \rightarrow e e \overline{\mu }}({{\textbf {p}}};\Delta t) = \sum _{{{\textbf {k}}} \ne {{\textbf {p}}},s_2,s_3} |\tilde{\mathcal {A}}^{ s_2 s_3}_{e \rightarrow e e\overline{\mu }}({{\textbf {k}}};t_i,t_f)|^2 \nonumber \\{} & {} + \sum _{s_2,s_3} |\tilde{\mathcal {A}}^{ s_2 s_3}_{e \rightarrow e e\overline{\mu }}({{\textbf {p}}};t_i,t_f)|^2 -\sum _{s_3} |\tilde{\mathcal {A}}^{r s_3}_{e \rightarrow e e\overline{\mu }}({{\textbf {p}}};t_i,t_f)|^2 \nonumber \\{} & {} = \sum _{{{\textbf {k}}} \ne {{\textbf {p}}},s_2,s_3} |\tilde{\mathcal {A}}^{ s_2 s_3}_{e \rightarrow e e\overline{\mu }}({{\textbf {k}}};t_i,t_f)|^2+\sum _{s_3} |\tilde{\mathcal {A}}^{r s_3}_{e \rightarrow e e\overline{\mu }}({{\textbf {p}}};t_i,t_f)|^2 \,.\nonumber \\ \end{aligned}$$
(95)

In other words, because of the Pauli principle, the vacuum should not carry the contribution with \({\varvec{k}}={\varvec{p}}\). Then, in the large-V limit

$$\begin{aligned} \mathcal {P}_{e \rightarrow e e \overline{\mu }}({{\textbf {p}}};\Delta t) \ {}= & {} \ V \sum _{s_2,s_3} \, \int \!\! \frac{\textrm{d}^3 {{\textbf {k}}}}{(2 \pi )^3} \, |\tilde{\mathcal {A}}^{ s_2 s_3}_{e \rightarrow e e \overline{\mu }}({{\textbf {k}}};t_i,t_f)|^2\nonumber \\{} & {} +\sum _{s_3} |\tilde{\mathcal {A}}^{r s_3}_{e \rightarrow e e\overline{\mu }}({{\textbf {p}}};t_i,t_f)|^2 \,. \end{aligned}$$
(96)

The first piece diverges and must be subtracted, while the second piece gives a finite contribution. Explicitly

$$\begin{aligned} \mathcal {P}_{e \rightarrow e e \overline{\mu }}({{\textbf {p}}};\Delta t)= & {} \frac{4 m_{e \mu }^2 Y_{{\textbf {p}}}^2 }{\left( \omega _{{{\textbf {p}}},e}+\omega _{{{\textbf {p}}},\mu }\right) ^2}\times \nonumber \\ {}{} & {} \times \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }+\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) , \end{aligned}$$
(97)

where

$$\begin{aligned} Y_{{{\textbf {p}}}}^2 \ = \ \sum _{{{\textbf {s}}}} \left( Y^{rs}_{{{\textbf {p}}}}\right) ^* Y^{rs}_{{{\textbf {p}}}} \,, \end{aligned}$$
(98)

and

$$\begin{aligned} Y_{{{\textbf {p}}}} \ = \ \frac{|{{\textbf {p}}}|}{\sqrt{4 \omega _{{{\textbf {p}}},e}\omega _{{{\textbf {p}}},\mu }}} \left( \sqrt{\frac{\omega _{{{\textbf {p}}},\mu }+m_\mu }{\omega _{{{\textbf {p}}},e}+m_e}}+\sqrt{\frac{\omega _{{{\textbf {p}}},e}+m_e}{\omega _{{{\textbf {p}}},\mu }+m_\mu }}\right) \,.\nonumber \\ \end{aligned}$$
(99)

Therefore, the total decay probability of \(\nu _e\) is

$$\begin{aligned}{} & {} \mathcal {P}^{e}_{D}({{\textbf {p}}}; \Delta t) = 4 m^2_{e \mu } \nonumber \\ {}{} & {} \quad \left[ \frac{W_{{{\textbf {p}}}}^2}{\left( \omega _{{{\textbf {p}}},e}-\omega _{{{\textbf {p}}},\mu }\right) ^2} \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) \right. \nonumber \\{} & {} \left. + \frac{Y_{{\textbf {p}}}^2 }{\left( \omega _{{{\textbf {p}}},e}+\omega _{{{\textbf {p}}},\mu }\right) ^2} \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }+\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) \right] \,. \end{aligned}$$
(100)

Note that all such decays (flavor transitions) are forbidden when \(t_i \rightarrow -\infty \), \(t_f \rightarrow +\infty \), unless \(m_e = m_\mu \). In particular, the last two processes are always forbidden for infinite time-intervals because of energy conservation. In fact, in that case, the three-dimensional \(\delta \)s would be substituted by deltas on the four-momentum, which strictly employ energy conservation. As we commented above, this is expected in analogy of what happens for unstable particles.

If we now define, following the notation of Ref. [15]

$$\begin{aligned} |U_{{\textbf {p}}}|= & {} W_{{{\textbf {p}}}}\frac{m_\mu -m_e}{\omega _{{{\textbf {p}}},e}-\omega _{{{\textbf {p}}},\mu }} \nonumber \\= & {} \sqrt{\frac{\left( \omega _{{{\textbf {p}}},e}+m_e\right) \left( \omega _{{{\textbf {p}}},\mu }+m_\mu \right) }{4\omega _{{{\textbf {p}}},e}\omega _{{{\textbf {p}}},\mu }}}\times \nonumber \\{} & {} \times \left( 1+\frac{|{{\textbf {p}}}|^2}{\left( \omega _{{{\textbf {p}}},e}+m_e\right) \left( \omega _{{{\textbf {p}}},\mu }+m_\mu \right) }\right) \,, \end{aligned}$$
(101)
$$\begin{aligned} |V_{{\textbf {p}}}|= & {} Y_{{{\textbf {p}}}}\frac{m_\mu -m_e}{\omega _{{{\textbf {p}}},e}+\omega _{{{\textbf {p}}},\mu }} \nonumber \\= & {} \ \sqrt{\frac{\left( \omega _{{{\textbf {p}}},e}+m_e\right) \left( \omega _{{{\textbf {p}}},\mu }+m_\mu \right) }{4\omega _{{{\textbf {p}}},e}\omega _{{{\textbf {p}}},\mu }}}\times \nonumber \\ {}{} & {} \times \left( \frac{|{{\textbf {p}}}|}{\omega _{{{\textbf {p}}},e}+m_e}-\frac{|{{\textbf {p}}}|}{\omega _{{{\textbf {p}}},\mu }+m_\mu }\right) \,, \end{aligned}$$
(102)

we can rewrite the probability in the form

$$\begin{aligned} \mathcal {P}^{e}_{D}({{\textbf {p}}}; \Delta t)= & {} \sin ^2 2 \theta \left[ |U_{{\textbf {p}}}|^2 \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) \right. \nonumber \\{} & {} \left. + |V_{{\textbf {p}}}|^2 \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }+\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) \right] . \end{aligned}$$
(103)

with \(\theta =m_{e \mu }/(m_\mu -m_e) \approx \sin \theta \). In the approximation we used this coincides with the oscillation probability (A11), firstly derived in [26]. This is a remarkable result because the method we adopted is quite different from the approach of Ref. [26] (which is briefly reviewed in appendix), and it does not rely on the construction of a flavor vacuum (see Eq. (A6)).

The second term on the r.h.s. of the Eq. (103) is the main correction with respect to the usual Pontecorvo oscillation formula (A16). The effect of such term is negligible for relativistic neutrinos, e.g. when \(m_\sigma /|{{\textbf {p}}}| \rightarrow 0\), while is maximal when \(|{{\textbf {p}}}| = \sqrt{m_e m_\mu }\) which, in our approximation is equivalent to \(|{{\textbf {p}}}| = \sqrt{m_1 m_2}\), i.e. when the momentum is of the order of neutrino masses. In such a case

(104)

where

(105)

Possible phenomenological implications in this regime could be studied in a cosmological context, how it will be discussed in the conclusions.

Let us now compute the survival probability \(\mathcal {P}^e_{S}({{\textbf {k}}},\Delta t)\). The zeroth order contribution gives \(\mathcal {P}^e_{S}({{\textbf {k}}},\Delta t)=1\). To find a non-trivial contribution to the survival probability we should compute the second order terms. In the present case it is useful to write

$$\begin{aligned}{} & {} U(t_i,t_f) = 1\hspace{-1mm}\textrm{I}-i \, m_{e \mu } \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x \,: \overline{\nu }_e(x) \nu _\mu (x)+\overline{\nu }_\mu (x) \nu _e(x): \nonumber \\{} & {} \quad - \ \frac{m^2_{e \mu }}{2} \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x_1 \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x_2 \, \, \mathcal {T}\Big [ \left( :\overline{\nu }_e(x_1) \nu _\mu (x_1) +\overline{\nu }_\mu (x_1) \nu _e(x_1):\right) \nonumber \\{} & {} \quad \left( : \overline{\nu }_e(x_2) \nu _\mu (x_2)+\overline{\nu }_\mu (x_2) \nu _e(x_2):\right) \Big ] \, + \,\ldots \,.\nonumber \\ \end{aligned}$$
(106)

The second order piece can be further expanded using Wick theorem:

$$\begin{aligned} U^{(2)}(t_i,t_f)= & {} -\frac{m^2_{e \mu }}{2} \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x_1 \nonumber \\{} & {} \times \int ^{t_f}_{t_i} \!\! \textrm{d}^4 x_2 \, \, \Big [:\overline{\nu }_e(x_1) \nu _\mu (x_1) \overline{\nu }_e(x_2) \nu _\mu (x_2):\ \nonumber \\{} & {} +:\overline{\nu }_e(x_1) \nu _\mu (x_1)\overline{\nu }_\mu (x_2) \nu _e(x_2): \nonumber \\{} & {} +:\overline{\nu }_\mu (x_1) \nu _e(x_1)\overline{\nu }_e(x_2) \nu _\mu (x_2):\ \nonumber \\{} & {} \ +:\overline{\nu }_\mu (x_1) \nu _e(x_1)\overline{\nu }_\mu (x_2) \nu _e(x_2): \nonumber \\{} & {} + 2 \, i \left( S^e_{\alpha \beta }(x_2-x_1) \,:\overline{\nu }^\beta _\mu (x_2) \nu ^\alpha _\mu (x_1):\right. \nonumber \\{} & {} \left. + S^\mu _{\alpha \beta }(x_2-x_1) \,:\overline{\nu }^\beta _e(x_2) \nu ^\alpha _e(x_1):\right) \Big ] \,, \end{aligned}$$
(107)

where

(108)

is the Dirac propagator. The employment of Wick theorem makes evident that the divergent vacuum contributions come from the terms containing the propagator, which can be neglected.

The survival process is

$$\begin{aligned} |\nu ^r_{{{\textbf {p}}},e}\rangle \ \rightarrow \ |\nu ^s_{{{\textbf {k}}},e}\rangle \,. \end{aligned}$$
(109)

Saving only up to linear terms in \(m_{e \mu }\), the amplitude can be written as

$$\begin{aligned} \mathcal {A}^{rs}_{e \rightarrow e}({{\textbf {p}}},{{\textbf {k}}}; t_i,t_f)= & {} \delta _{{{\textbf {k}}}, {{\textbf {p}}}} \delta _{rs}+\frac{1}{2}\mathcal {A}^{(2) r s}_{e \rightarrow e}({{\textbf {p}}},{{\textbf {k}}}; t_i,t_f) \,,\nonumber \\ \end{aligned}$$
(110)

where \(A^{(2) r s}_{e \rightarrow e}({{\textbf {k}}},{{\textbf {p}}};t_i,t_f)\) is the second-order piece, which goes as \(m_{e \mu }^2\). Taking the square, and summing over the final momenta and helicities, the only pieces which cannot be disregarded are the one which go as \(m_{e \mu }^2\) or lower powers, i.e.

$$\begin{aligned} \mathcal {P}^{e}_{S}({{\textbf {p}}}; \Delta t) \ {}= & {} \ \sum _{{{\textbf {k}}},s}\mathcal {A}^{rs}_{e \rightarrow e}({{\textbf {p}}},{{\textbf {k}}}, t_i,t_f) \nonumber \\\approx & {} \ 1+ 2 \, \Re e \left( \tilde{\mathcal {A}}^{(2)}_{e \rightarrow e}({{\textbf {p}}};t_i,t_f)\right) \,, \end{aligned}$$
(111)

with

$$\begin{aligned} \tilde{\mathcal {A}}^{(2)}_{e \rightarrow e}({{\textbf {p}}};t_i,t_f) \ \equiv \ \sum _{{{\textbf {k}}},s} A^{(2) r s}_{e \rightarrow e}({{\textbf {p}}},{{\textbf {k}}};t_i,t_f) \,. \end{aligned}$$
(112)

Explicitly one finds

$$\begin{aligned} \mathcal {P}^{e}_{S}({{\textbf {p}}}; \Delta t)= & {} 1-\sin ^2 2 \theta \nonumber \\ {}{} & {} \times \left[ |U_{{\textbf {p}}}|^2 \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }-\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) \right. \nonumber \\{} & {} \left. + |V_{{\textbf {p}}}|^2 \sin ^2\left( \frac{\left( \omega _{{{\textbf {p}}},\mu }+\omega _{{{\textbf {p}}},e}\right) \Delta t}{2}\right) \right] \,. \end{aligned}$$
(113)

Then

$$\begin{aligned} \mathcal {P}^{e}_{D}({{\textbf {p}}}; \Delta t) + \mathcal {P}^{e}_{S}({{\textbf {p}}}; \Delta t)\ = \ 1 \,, \end{aligned}$$
(114)

as expected.

6 Conclusions

In this paper we have shown that neutrino oscillations can actually be described as a standard perturbative QFT, with the mixing handled by the interaction picture. A 0+1D toy model, “scalar” neutrinos and their “realistic” Dirac fermion counterparts can all be described by this approach, the various transition amplitudes can be related via Feynman diagrams and the differences reduced to spin-statistics differences. In the fermionic case, the same oscillation formula as in non-perturbative flavor-Fock space approach [26], is here independently recovered, within the approximations involved in the perturbative calculation. In particular, we find a term which depends on the sum of the frequencies, in addition to the usual Pontecorvo term which only involves their difference. Moreover, the pre-factors in front of such terms are exactly the Bogoliubov coefficients derived in Ref. [15].

This description sheds new light on the long-standing arguments in the literature as to the “true basis” of neutrino Fock space. As in other quantum field theories, provided interacting and free states have the same quantum numbers, the Kallen-Lehmann representation can be used to convert between the two perturbatively. In such a picture, experimental data selects the physical degrees of freedom, to be related to the Lagrangian ones via Renormalization. Since interaction happens via weak charge, the perturbation series defined in this work shows that neutrino physical states can be consistently defined.

The methods developed in this work can be applied to a calculation of observables beyond tree level which systematically includes oscillation effects. An obvious candidate is the magnetic dipole moment of the neutrino, which so far has only been obtained by effective field theory techniques [43]. A proper calculation will include the interaction included here (which corrects the magnetic mass) with the same footing as the W-lepton bubble (which provides the electric current generating the magnetic moment). Given the non-trivial momentum dependence of QFT neutrino oscillation effects noticed in [15] the interplay of all these effects could yield surprises, which might be observable [44].

There are however subtleties not present in other usual quantum field theories. Unlike other interactions, neutrino oscillations are 2-field operators, analogous to coupled harmonic oscillators. The vacuum of the theory should therefore be analytically obtainable, and it is well known that perturbative terms can not capture such vacuum corrections. As a related issue, the self-interaction due to the vacuum self energy will result in a “mass” correction to the flavor state, tilting the perturbative series examined in this work. Since the theory examined here is perfectly Lorentz invariant, violations derived using the quantum field theory approach to fermion mixing (see e.g. [45]) can only be a feature of the vacuum and to properly assess such effects a perturbative series must be derived from such a vacuum. This is left to a forthcoming work.

Moreover, it is interesting to stress that the physical corrections to the usual neutrino oscillation formula (A16) could be relevant in the measurements of the so-called cosmic neutrino background (CNB). This will be studied, e.g., by the PTOLEMY experiment [46,47,48,49]. The basic idea of such experiment is to detect the relic neutrinos of the CNB (which decoupled around one second after the Big Bang) through the neutrino capture by tritium \(\nu _e +{}^3 \textrm{H}\rightarrow e^- + {}^3\textrm{He}\). With the standard approach to neutrino oscillations, the rate of the capture process reads [49]

$$\begin{aligned} \Gamma _{CNB} \ = \ \sum _j \, |U_{e j}|^2 \, \bar{\sigma } \, v_\nu \, f_{e,i} \, n_0 \,, \end{aligned}$$
(115)

where U is the mixing matrix, \(v_\nu \) is the neutrino velocity as measured at Earth, \(\bar{\sigma }\) is the average cross-section, \(n_0\) is the average neutrino number density on large scale and \(f_{e,i}\) are the clustering factors. We expect that the analysis of this paper modifies the above expression, with the inclusion of the coefficients \(|U_{{\textbf {p}}}|\) and \(|V_{{\textbf {p}}}|\). This could be easily understood looking at Ref. [21], where the rate of \(\beta \)-decay \({}^3 \textrm{H} \rightarrow e^- + {}^3\textrm{He}+\nu _e\) was computed within the flavor Fock-space approach. As discussed in the text (see Eqs. (104)–(105)), the corrections in that case, and then in the present case, become relevant in the non-relativistic regime. However, a quantitative analysis requires to consider the three-flavor phenomenology and it will be postponed to a forthcoming work.

Finally, another phenomenon which is suited to be investigated within the approach here introduced are chiral oscillations [50,51,52], which again are relevant in the non-relativistic regime and have been recently discussed in connection with CNB [53]. Work is in progress in this direction.