Neutrinos in curved spacetimes: particle mixing and flavor oscillations

We present a quantum field theoretical approach to the vacuum neutrino oscillations in curved space, we analyze the non--trivial interplay between quantum field mixing and field quantization in curved space and derive new oscillation formulae. We compute the formulae explicitly in the spatially flat FLRW metrics for universes dominated by a cosmological constant and by radiation. We evaluate the transition probabilities in the Schwarzschild black hole metric, and we show that the Hawking radiation affects the oscillations of neutrinos. We show that our results are consistent with those of previous analyses when the quantum mechanical limit is considered.

We present a quantum field theoretical approach to the vacuum neutrino oscillations in curved space, we analyze the non-trivial interplay between quantum field mixing and field quantization in curved space and derive new oscillation formulae. We compute the formulae explicitly in the spatially flat FLRW metrics for universes dominated by a cosmological constant and by radiation. We evaluate the transition probabilities in the Schwarzschild black hole metric, and we show that the Hawking radiation affects the oscillations of neutrinos. We show that our results are consistent with those of previous analyses when the quantum mechanical limit is considered.

I. INTRODUCTION
Since they were theoretically proposed by Pauli [1], neutrinos have proven to be among the most enigmatic particles in the universe. Until the discovery of flavor oscillations [2,3], whose theory was pioneered by Pontecorvo [4,5], neutrinos were believed to be massless. Today it is accepted that neutrinos are massive particles, and that they oscillate among three flavors ν e , ν µ , ν τ corresponding to the companion charged leptons e, µ, τ . This peculiarity renders neutrinos unique among the known elementary particles and puts them beyond the scope of the standard model of particles [6]. In many respects, neutrinos are forerunners of a new physics, as several issues, including the origin of their mass [7] and their fundamental nature [8], are still open to the present day.
On the other hand, the relevance of neutrinos in astrophyisical and cosmological contexts has grown dramatically during the last years. They figure as a valuable source of information, along with gravitational waves and electromagnetic radiation, in the ever-growing field of multi-messenger astronomy [9]. The study of neutrinos of astrophysical origin can indeed provide fundamental insights on the source that produced them. In addition, neutrinos are expected to play an important role in the first phases of the universe [10,11], and the detection of the cosmic neutrino background, pursued in experiments as PTOLEMY [12], could represent an essential test for the standard cosmological model [11]. Mass varying neutrinos have also been proposed as a possible explanation for Dark Energy [13].
This state of affairs requires a careful investigation of neutrino oscillations on a curved spacetime. The topic has been discussed in several works , where it was found that gravitational fields may alter both the oscillations in vacuum and in matter [14][15][16].
Here we wish to go beyond the heuristic treatment of ref. [16], and present a quantum field theoretical approach, based on the field quantization in curved space-time, to evaluate the effects of gravitational fields on neutrino oscillations. We derive general oscillation formulae for flavor fields in curved space-time, which represent our main result. We discuss the particle interpretation of the fields in presence of gravity and study how the mixing changes when moving from a mass field representation to another. We demonstrate the invariance of local observables, which are represented by expectation values on flavor states of local operators constructed from the flavor fields. We show that the oscillation probabilities, on the other hand, do in general depend on the representation of the mass fields, since they are not a local observable and involve the comparison between particles in different spacetime regions. We establish the conditions which have to be satisfied in order that the resulting transition probabilities are invariant under changes of mass field representation.
We also compute explicitly the oscillation formulae for two examples of spatially flat Friedmann-Lemaitre-Robertson-Walker spacetimes, corresponding to a cosmological constant-dominated and a radiation-dominated universe respectively. In these cases, exact analytical solutions to the Dirac equation are available, and the formalism here introduced can be applied directly. Moreover, we give an estimation of the oscillation formulae for neutrinos propagating from the past infinity to the future infinity in a stationary Schwarzschild spacetime. We introduce a method to extract the oscillation formulae on spacetimes with asymptotically flat regions without resorting to the exact solutions of the Dirac equation. We then employ this strategy to compute the formulae on the Schwarzschild black hole spacetime, for neutrinos propagating from the past infinity to the future infinity. We show how the Hawking radiation is naturally embedded in the resulting transition probabilities.
Our results generalize those of the previous treatments [16], and are consistent with the latter when the suitable limits are considered. In our computation, for simplicity, we limit our analysis to the vacuum oscillations, therefore considering the sole effect of gravity.
The paper is organized as follows: in section II we provide the setting for the description of the mass fields in curved space; in section III we develop field mixing and find the oscillation probabilities in curved spacetime, with a thorough analysis of their features; in section IV we apply the formalism to some spacetimes of interest, including the spatially flat FLRW metric for a radiation-dominated universe and for a cosmological constant-dominated universe, and the Schwarzschild black hole metric, where we show the impact of the Hawking effect on neutrino oscillations; finally in section V we draw our conclusions.

II. MASS NEUTRINO FIELDS IN CURVED SPACE
To evaluate the oscillation formulae for neutrinos on a curved spacetime, it is necessary to consider both the effects of curvature and mixing on the (free) mass fields. Let M be a globally hyperbolic spacetime, and let τ ∈ R label a foliation of M by Cauchy surfaces. Consider the tetrad fields e µ a (x) satisfying η ab e µ a (x)e ν b (x) = g µν (x). Here η ab ≡ diag(1, −1, −1, −1) is the Minkowski metric tensor, while g µν (x) is the contravariant metric g µν (x)g νρ (x) = δ µ ρ on M in a given coordinate system. The massive neutrino fields satisfy the Dirac equations: where γ µ (x) = e µ a (x)γ a , γ a being the usual flat space Dirac matrices, and D µ = ∂ µ − i 4 ω ab µ σ ab . The spin connection is defined as ω ab µ = e a ν Γ ν ρµ e ρb + e a ν ∂ µ e νb , whereas σ ab are the commutators of flat Dirac matrices σ ab = i 2 [γ a , γ b ]. In equation (1), the index i = 1, 2, ..., N ranges over the number of neutrino species N . For the sake of simplicity we focus on the case N = 2, though the generalization to N = 3 is straightforward. In general equation (1) cannot be solved exactly. Even if one is able to find exact solutions, these do not play the same prominent role as their flat spacetime counterpart. It is well-known, indeed, that the positive frequency solutions of equation cannot be defined univocally, and that, consequently, there is no natural (nor unique) particle interpretation for the corresponding Quantum Field Theory [17,18]. Nevertheless, the canonical quantization of the Dirac field proceeds along the same lines as in Minkowski spacetime.
To perform a field expansion, one must find a set of positive ζ k,i and negative ξ k,i frequency solutions for each of the equations (1). In general the bipartition of the solutions to eq. (1) makes sense only locally, while there is no natural global definition of positive and negative frequency modes. Anyway, one is free to choose a set of modes [19] {ζ k,i , ξ k,i } , deemed to be positive/negative frequency modes according to some specified observer, and expand the field with respect to them, provided that they form a complete (and orthonormal) set of solutions under the inner product with a i , b i any solution to equation (1) with mass m i andb i = b † i γ 0 (x). Here dΣ µ (τ ) = n µ (τ )dV τ denotes the volume element on the surface τ with unit timelike normal n µ (τ ). This has to hold separately for each i = 1, 2. As it is easy to prove, for a i , b i solutions of the (same) Dirac equation, the inner product (2) does not depend on the hypersurface chosen for the integration. In particular, it is independent on the foliation by Cauchy hypersurfaces employed. The fields can then be expanded as with the operator coefficients γ k,s,i , ǫ k,s,i satisfying the usual canonical anticommutation relations, k momentum index and s helicity index. The annihilators are also required to anticommute for i = j, and, in particular {γ k,s,i , γ † k,s,j } = δ ij , {ǫ k,s,i , ǫ † k,s,j } = δ ij . In equation (3) we prefer to keep any space-time dependence within the modes, for ease of treatment with a general metric. The expansions (3) define the mass Hilbert space H m = H 1 ⊗ H 2 , which is constructed out of the vacuum |0 m = |0 1 ⊗ |0 2 . Here |0 i is defined, as usual, by γ k,s,i |0 i = 0 = ǫ k,s,i |0 i for each k, s, i.
As hinted above, the field expansions (3) are somewhat arbitrary, as opposed to the flat spacetime case, where there is no ambiguity in the definition of positive and negative frequency modes. Any other basis {ζ k,s,i ,ξ k,s,i } can be used to expand the fields ψ i = k,s (γ k,s,iζk,s,i (x) +ǫ † k,s,iξ k,s,i (x)). Since both the sets {ζ k,s,i , ξ k,s,i } and {ζ k,s,i ,ξ k,s,i } form a basis for the space of solutions of eqs. (1), one can write the modes of a set in terms of the other, for each i: where Γ k ′ ,s ′ ;ks;i = (ζ k ′ ,s ′ ,i , ζ k,s,i ) = (ξ k,s,i ,ξ k ′ ,s ′ ,i ) and Σ k ′ ,s ′ ;ks;i = (ζ k ′ ,s ′ ,i , ξ k,s,i ) = −(ζ k,s,i ,ξ k ′ s ′ ,i ). This is a fermionic Bogoliubov transformation, for which q,r Γ * k,s;q,r;i Γ k ′ ,s ′ ;q,r;i + Σ * k,s;q,r;i Σ k ′ ,s ′ ;q,r;i = δ k,k ′ δ s,s ′ for each i. The corresponding relation between the two sets of annihilators is given bỹ It is often the case that the Bogoliubov coefficients Γ k,s;k ′ ,s ′ ;i , Σ k,s;k ′ ,s ′ ;i can be written as Γ k,s;k ′ ,s ′ ;i = δ k,k ′ δ s,s ′ Γ k,i , Σ k,s;k ′ ,s ′ ;i = δ k,k ′ δ s,s ′ Σ k,i , with Γ k,i and Σ k,i depending on k alone. In this occurrence, they admit the parametrization Γ k,i = e iη k,i cos(θ k,i ) ,Σ k,i = e iφ k,i sin(θ k,i ), with η k,i , φ k,i , θ k,i real functions of k. We remark that the Bogoliubov transformations (5) can be recast in terms of the generators The maps J i :H i → H i interpolate between the Fock spaces H i built from the γ k,s,i , ǫ k,s,i and the Fock spaceH i built from theγ k,s,i ,ǫ k,s,i . In particular, one has for the vacuum states |0 i = J −1 i |0 i . As for the untilded representation, the mass Hilbert space in the tilded representation is the tensor productH m =H 1 ⊗H 2 . It is convenient to define a unique generator of Bogoliubov transformations J :H m −→ H m onH m as the tensor product J = J 1 ⊗ J 2 . Theñ for i = 1, 2. The expansions of the two fields ψ 1 and ψ 2 must be compatible with each other, i.e., each of the modes ζ k,s,2 , ξ k,s,2 must be obtainable from the corresponding modes ζ k,s,1 , ξ k,s,1 by the substitution m 1 ↔ m 2 , and vice-versa. In the context of mixing, this ensures that the same kind of particle, described by the same set of quantum numbers, is being mixed. This, of course, does not undermine the arbitrariness in the choice of the modes; these can be any complete set of solutions to the Dirac equation, provided that the same choice is made for the two fields.

III. NEUTRINO MIXING AND OSCILLATION FORMULAE IN CURVED SPACE-TIME
In this Section, we show new oscillation formulae for flavor fields in curved space-time and we present general considerations on the infinitely many unitarily inequivalent representations of the canonical anticommutation relations which characterize the quantization of mixed fields and of fields in curved space.

A. Oscillation Formulae
As discussed above, the QFT of free Dirac fields in curved space is characterized by infinitely many unitarily inequivalent representations of the canonical anticommutation relations. The phenomenon of mixing, even in Minkowski space, suffers from an analogous ambiguity, in that the flavor and the mass representations are unitarily inequivalent [20]. The effects of such inequivalence have been analyzed in flat space time [21] and the possibility to reveal them in experimental setup has been recently proposed [22]. Let us start by fixing the mass field expansions (3) and describe the mixing in a given representation of the mass fields. The flavor fields are defined as ψ e = cos(θ)ψ 1 + sin(θ)ψ 2 and ψ µ = cos(θ)ψ 2 − sin(θ)ψ 1 with θ the (2-flavor) mixing angle. Just like the Bogoliubov transformations (7), the rotation to flavor fields can be cast in terms of a generator I θ (τ ). This is given by where the scalar products (ψ i , ψ j ) do depend on the hypersurface chosen for the integration, since they are solutions to different Dirac equations. Then, by definition, the flavor fields are expressed as If we let the generator (8) act on the mass annihilators, we obtain the flavor annihilators for curved space And similar for γ k,s,µ (τ ), ǫ k,s,e (τ ), ǫ k,s,µ (τ ). The Bogoliubov coefficients are provided by the inner products of the solutions to the curved space Dirac equation with mass m 1 and m 2 , that is, Λ q,r;k,s (τ ) = (ζ q,r,2 , ζ k,s,1 ) τ = (ξ k,s,1 , ξ q,r,2 ) τ and Ξ q,r;k,s (τ ) = (ζ k,s,1 , ξ q,r,2 ) τ = −(ζ q,r,2 , ξ k,s,1 ) τ . The mixing coefficients always satisfy Since the mass expansions are compatible, the mixing coefficients are often diagonal, namely of the form with Λ k,s (τ ), Ξ k,s (τ ) depending on k and s alone [23]. Exceptions to this arise when we consider expansions of the mass fields in terms of modes labelled by the energy. In such a case, the mixing coefficients are non-diagonal and different from zero for each k, s, τ , and respectively. The mass and the flavor representations are unitarily inequivalent. For each τ one has a distinct flavor Fock space In order to define the transition probabilities, we observe that the total Lagrangian is invariant under U (1) gauge transformations. Therefore the total charge It is then meaningful to define the transition probabilities as Here ρ, σ = e, µ, the state |ν ρ,k,s (τ 0 ) = γ † k,s,ρ (τ 0 ) |0 f (τ 0 ) is the state with a single neutrino of flavor ρ, momentum k and helicity s on the reference hypersurface τ = τ 0 . The second term on the rhs of (15) is just the implementation of the normal ordering with respect to |0 f (τ 0 ) . By construction P e→e k,s (τ ) + P e→µ k,s (τ ) = 1 and P µ→e k,s (τ ) + P µ→µ k,s (τ ) = 1 for each τ . A straightforward calculation yields, in the general case (accounting for both diagonal and non-diagonal mixing coefficients) the result Equation (16) is the central result of the paper. When equations (12) hold, this reduces to In both cases one has B. Mixing on a curved background and gravity-induced ambiguity in the particle interpretation Up to now, we have worked within a fixed, but arbitrary, representation of the mass fields. The question arises about the other possible representations, and how the mixing changes when moving from a representation to another. For the definition (15) to make sense, we must determine if and how the probabilities vary when the mass representation is changed. We take as a guideline the principle of covariance, so that the local physical observables should be independent of the underlying representation. In moving from a given representation {γ 1 , ǫ 1 }, {γ 2 , ǫ 2 } to another {γ 1 ,ǫ 1 }, {γ 2 ,ǫ 2 }, we know how to connect the mass Fock spaces, namely via the generator (7) J −1 : H m →H m . For each mass representation, we can proceed as we did above and build the corresponding flavor annihilators and flavor spaces H f (τ ),H f (τ ), together with the mixing generators I θ (τ ) : It is useful, at this point, to determine the relations among the mixing coefficients Λ(τ ), Ξ(τ ) and Λ(τ ),Ξ(τ ) that appear in the explicit form of the two generators I θ (τ ) andĨ θ (τ ). By definition we havẽ Here the first equality is just the definition ofΛ(τ ), the second follows from the Bogoliubov transformations (5). By using the properties of the inner product (2), in the general case (again, accounting for both the diagonal and non-diagonal mixing coefficients), we obtaiñ and, finally, from the definition of Λ(τ ) and Ξ(τ ), we havẽ Similarly we haveΞ When equations (12) hold for both the representations {q, r}, {q ′ , r ′ }, the equations reduce tõ The equations (21,22,23,24) provide an explicit relation between I θ (τ ) andĨ θ (τ ), and show how the mixing coefficients change, in moving from a mass representation to another, in order to ensure covariance. In particular, the tilded coefficients turn out to be a linear combination of the untilded coefficients weighted by the coefficients of the Bogoliubov transformations between the two mass representations. A slightly modified version of the eqs. (21) and (22) will be expedient in the calculation of the transition probabilities in a number of interesting cases.
It remains to establish how the flavor operators γ ρ (τ ), ǫ ρ (τ ) and the flavor vacuum |0 f (τ ) transform under a change of mass representation. We focus on the vacuum state |0 f (τ ) ∈ H f (τ ). First we employ the generator I θ (τ ) : H f (τ ) → H m to get the mass vacuum |0 m . Then we apply the generator of Bogoliubov transformations (7) J −1 : H m →H m to obtain |0 m . Finally, the generator of mixing in the tilde representationĨ −1 θ (τ ) :H m →H f (τ ) is employed to get |0 f (τ ) . We conclude that the two flavor vacua are related by the transformation .

C. Transition probabilities and the mass representation
The oscillation probabilities are not a local observable, since they involve the comparison between particles at different values of τ , as it is evident from the definition (15). In general these quantities do depend on the representation of the mass fields, and this is because distinct representations might assign a different meaning to the quantum numbers k, s. For example, one might consider two expansions of the mass fields, one in terms of plane waves, labelled by the three-momenta {k k k}, and one in terms of localized wave packets, labelled by a suitable set of quantum numbers {q}. It is clear that the two expansions describe particles with different physical properties; the first describes particles with definite momentum, the second describes particles for which momentum and position are definite to some extent. Therefore the probabilities P ρ→σ k,s and P ρ→σ q,s refer to the oscillations of different particles, and have a different interpretation. It would make no sense, in such a case, to require the equivalence of the two. This, of course, would be true even in flat space.
In the most general case, as the quantum numbers k, s and k ′ , s ′ have a different phyiscal meaning, the probabilities P ρ→σ k,s andP ρ→σ k ′ ,s ′ have different interpretations. Different representations of the mass fields do indeed assign a different meaning to such indices, so that the probabilities (15) have no invariant meaning. In order to make sense of the probabilities in eqs. (15) in the most general case, a representation of the mass fields must be fixed on the grounds of physical relevance. When the underlying spacetime M possesses non trivial symmetries, as time translational invariance or spherical symmetry, there is no doubt that the representation should be fixed so to take them into account. In these cases "good quantum numbers" are suggested by the symmetries themselves (for instance, the energy ω for stationary metrics, the angular momentum l, m for spherically symmetric spacetimes). In any case, the probabilities in a given mass representation can always be related to the probabilites in any other mass representation with the aid of equations (23) and (24). As a final remark, we stress that the issue discussed here has nothing to do with the diffeomorphism covariance of the theory. All the probabilities (17) are (generally covariant) scalars, as it is evident from the definitions.

IV. NEUTRINO OSCILLATION FORMULAE IN FLRW METRICS AND IN PRESENCE OF A SCHWARZSCHILD BLACK HOLE
In this section we apply the formalism developed above to some cases of interest. After an analysis of the flat space limit, we consider two cosmologically relevant FLRW metrics, corresponding to a cosmological constant-dominated and a radiationdominated universe respectively. In these cases, exact analytical solutions to the Dirac equation are available, and it is possible to employ equation (16) directly. We then introduce a method to extract the oscillation formulae on spacetimes with asymptotically flat regions without resorting to the exact solutions of the Dirac equation. We employ this strategy to compute the formulae on the Schwarzschild black hole spacetime, for neutrinos propagating from the past infinity to the future infinity. We show how the Hawking radiation is naturally embedded in the resulting transition probabilities.

A. Flat spacetime limit
As a first, trivial, application of the formulae (17), let us check the flat spacetime limit. We can see at once that the equations (17) reduce to the ordinary oscillation formulae. Indeed, in this case, we can choose the cauchy hypersurfaces to be the t = constant surfaces in a given Minkowskian coordinate system, while the modes {ζ k k k,s,i (x), ξ k k k,s,i (x)} are just the plane wave solutions to the flat Dirac equation with definite momentum, so that Λ q,r;k,s (t) → U q,r;k,s (t) and Ξ q,r;k,s (t) → V q,r;k,s (t), where U q,r;k,s and V q,r;k,s are the usual mixing coefficients in flat space [20]. Assuming, without loss of generality, k along the z direction, the helicity indices decouple U q,r;k,s = δ 3 (k k k − q q q)δ r,s U k k k , V q,r;k,s = δ 3 (k k k − q q q)δ r,s (−1) s V k k k . Since U k k k (t) = U k k k (0)e i(ω k,2 −ω k,1 )t and V k k k (t) = V k k k (0)e i(ω k,2 +ω k,1 )t , we get Substituting this result in eqs. (17) yields the flat space oscillation formulae, which further reduce to the Pontecorvo oscillation formulae in the quantum mechanical limit (V k k k = 0). Flat space also offers the possibility to illustrate some points discussed above in the simplest possible context. One might well expand the mass fields in terms of modes with definite energy and angular momentum ζ ω,κj ,mj ;i , ξ ω,κj ,mj ,s;i instead of considering modes with definite cartesian three-momentum k k k. The former shall be suitable combinations of spherical spinors [25]. An interesting aspect, is that in such a representation, the mixing coefficients are no longer diagonal. Indeed one has with ∆m 2 = m 2 2 − m 2 1 and and similar for Ξ, where the exponential is e i(ω ′ +ω)t and |U ω,ω ′ | is replaced by Here the quantum number κ j refers to a relativistic generalization of the spin-orbit operator, which enters the Dirac equation in spherical coordinates [25,26], and takes into account both the orbital and spin angular momentum. The index m j refers to one component of the total angular momentum J J J, and has not to be confused with the masses m i . Without delving into the details of the calculation, the result of equation (34) can be understood as follows. The modes, apart from a normalization constant, are given by where Ω κj ,mj (θ, φ) are spherical spinors, and P κj (λ i r) are radial functions of the product λ i r, with the radial momentum λ i = ω 2 − m 2 i . These are solutions to the radial part of the Dirac equation, which turns out to be a Riccati-Bessel equation [27]. The functions P κj (λ i r) are combinations of spherical bessel functions j n of the form P κj (λ i r) = rj κj (λ i r). In computing the inner products, the radial integration drP * κj ;2 (λ 2 r)P κj ;1 (λ 1 r) will produce a factor δ λ2,±λ1 , because of the closure relation satisfied by the spherical Bessel functions. Since λ 2 = ω ′ 2 − m 2 2 and λ 1 = ω 2 − m 2 1 , this will give rise to the delta factor appearing in (34). Notice that |U ω,ω ′ |, |V ω,ω ′ | are numerically the same as the usual flat space coefficients |U k k k |, |V k k k |, when This shows why the mixing coefficients Λ, Ξ are not generally diagonal, and the flexibility of the formalism we have employed. Indeed, the non-diagonal coefficients automatically ensure that the flavor operators γ ω,ρ , ǫ ω,ρ take into account the mass difference, involving operators with distinct energies ω, ω ′ for the fields ψ 1 and ψ 2 [28]. In flat space, the shift between the two representations is actually of no use. However, in a non-trivial framework, the versatility of the formalism is essential, as there are instances in which the cartesian components of the momentum k x , k y , k z are useless, while the "spherical" quantum numbers ω, l, m are well defined.

B. Expanding universe with exponential growth of the scale factor
The simplest non-trivial application is to spatially flat Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetimes. Consider the metric ds 2 = dt 2 − a 2 (t)(dx 2 + dy 2 + dz 2 ) with an exponential expansion a(t) = e Ht , H = constant. This is well-suited to describe a homogeneous, spatially flat and isotropic universe dominated by a cosmological constant. The normalized solutions to the Dirac equation for this metric were derived in [29]. Assuming, without loss of generality, the momentum k to be along the z direction, the helicities decouple as in flat space. Choosing the Cauchy surfaces as the surfaces with t = constant, or equivalently with a(t) = constant, the mixing coefficients read C. Expanding universe dominated by radiation Here we consider the FLRW metric for a radiation dominated universe a(t) = a 0 t 1 2 . Notice that since a(t) has to be adimensional, a 0 has dimension [t] − 1 2 = [m] 1 2 . As before, without loss of generality, we assume the neutrino momentum k along the z direction to decouple the helicities, and consider a foliation by the t = constant hypersurfaces. The solutions to the Dirac equation for this metric are again found in [29], and yield the mixing coefficients Λ k,s;q,r (t) = δs,rδ 3 where W κ,µ (z) are the Whittaker functions [27] and κ j = 1 4 1 + 2ik 2 a 2 0 mj for j = 1, 2. Insertion in eqs. (17) gives the transition probabilities (t 0 , t > 0) P e→µ k,s (t) = 2 cos 2 (θ) sin 2 (θ) 1 + ℜ

D. Spacetimes with asymptotically flat regions
The FLRW spacetimes considered above are among the few non-trivial metrics for which the Dirac equation (1) 17) is prohibitive. Nevertheless, if one is able to provide the relation between the two sets A, B, in the form of a Bogoliubov transformation, it is possible to derive oscillation formulae for the propagation from Ω A to Ω B . Assume that Ω A = τ ≤τA Σ τ and Ω B = τ ≥τB Σ τ for some values τ B > τ A , and consider τ 0 ≤ τ A , τ ≥ τ B . Suppose also that the B modes are given in terms of the A modes as Here the Bogoliubov coefficients Γ k ′ ,s ′ ;k,s;i , Σ k ′ ,s ′ ;k,s;i are again provided by the inner products (u A i , u B i ), (u A i , v B i ), yet their significance is slightly different from those in (4). While equation (4) describes a general and arbitrary change of basis, the transformation of equation (45) is dictated by the circumstance of having a natural choice for the modes in Ω A and Ω B , with a well-defined physical meaning. Indeed, the Bogoliubov coefficients take into account the effect of the curvature in the intermediate region. Then we can specialize equations (21), (22) to get For τ < τ A Λ A (τ ), Ξ A (τ ) are trivial, and vice-versa, for τ > τ B , Λ B (τ ), Ξ B (τ ) are trivial. Choosing, for instance, the B representation, one would have a trivial expression for Λ B (τ ), Ξ B (τ ) and one can make use of eqs. (47), (48) to obtain , which are also trivial. Then one can plug Λ B (τ, τ 0 ) and Ξ B (τ, τ 0 ) in equation (17) to obtain P e→µ k,s→q,r (τ ) for τ ≥ τ B and the reference hypersurface τ 0 ≤ τ A .

Scwharzschild black hole
As a realization of this scheme, consider the (static) Schwarzschild metric and two sequences of spacelike hypersurfaces [30] {Σ + n } n∈N , {Σ − n } n∈N approaching, respectively, the future I + and the past null infinity I − as n → ∞. We require that Σ + n and Σ − n be Cauchy surfaces respectively for the causal past J − (I + ) of I + and the causal future J + (I − )of I − for each n. For n large enough, these surfaces span an approximately flat portion of the Schwarzschild spacetime. On the surfaces Σ − n , as n approaches infinity, we expand the massive fields in terms of the incoming solutions, with frequency defined with respect to the schwarzschild time t and with definite angular momentum, These modes reduce to the flat space solutions (37) and (38) as I − is approached. Omitting the irrelevant angular and spin quantum numbers, as n → ∞, we get with |U ω,ω ′ | and |V ω,ω ′ | flat space spherical mixing coefficients as defined in (35) and (36), and ϕ − (ω, n), ρ − (ω, n) phase factors depending on ω and n. A similar reasoning can be carried on for the outgoing modes emerging at I + , ζ OUT ω,κj,mj ;i (t, r, θ, φ) , ξ OUT ω,κj ,mj;i (t, r, θ, φ), so to yield as n → ∞ It is understood that these equations hold as long as the spacetime is stationary (eternal black hole). For a body that collapses to a Schwarzschild black hole, as considered in the original paper by Hawking [32], we expect a slightly modified version of the Bogoliubov transformations, comprising non-diagonal thermal coefficients Γ ω,ω ′ ;i and Σ ω,ω ′ ;i with ω = ω ′ . We pick the ingoing representation to calculate the probabilities (of course, the outgoing representation yields the same result, as the Bogoliubov transformations are diagonal), and employ eqs. (47, 48) to obtain, with F H (ω) = and Choosing as reference hypersurface Σ − m for large m, we can now compute the probabilities (17) for a neutrino propagating from Σ − m to Σ + n , i.e. from I − to I + in the limit m, n → ∞. We find, for m, n → ∞ − F H (ω)F H (ω ′ ) |U ω;ω ′ | 2 cos(∆ + ω;m,n ) + |V ω;ω ′ | 2 cos(Φ + ω;m,n ) .

E. Quantum mechanical limit
It is a general feature of equations (17) that when all the quantum field theoretical effects are negligible, the oscillation formulae are modified only for a phase factor. Indeed, when one can neglect Ξ k,s in equation (17), one immediately has |Λ * k,s (τ )| = 1, (|Λ ω,ω ′ (τ )| = 1 respectively, in the non-diagonal case) for each τ from equation (12). Then the product Λ k,s (τ 0 ) * Λ k,s (τ ), (Λ ω,ω ′ (τ 0 ) * Λ ω,ω ′ (τ ) respectively) is just a phase e iϕ(τ0,τ ) and the net effect is a phase shift with respect to the Pontecorvo oscillation formulae, consistently with previous results obtained in a heuristic treatment [16]. Of course, the explicit value of the phase ϕ(τ 0 , τ ) depends on the metric and the surfaces τ considered, as well as on the mode expansion chosen for the mass fields. When the gravitational fields are weak enough, the phase can be computed by means of geometrical optics considerations [16,33].

V. CONCLUSIONS
We have developed a quantum field theoretical approach to the vacuum neutrino oscillations in curved space, discussing the transition probabilities, and their behaviour under changes of mass representation. We have analyzed the non-trivial interplay between quantum field mixing and field quantization in curved space, and have found that the former has a remarkably richer structure when compared to its flat space counterpart. In particular, the formalism has to be versatile enough to deal with the existence of infinite unitarily inequivalent representations of the mass fields, which is to say that no preferred notion of particle does generally exist. In the spirit of general covariance, we have determined the effect on flavor fields of a shift in the expansion of the mass fields, and established under which conditions the resulting transition probabilities are the same. We have then computed the oscillation formulae in three example metrics, including two FRLW spacetimes and the static Schwarzschild black hole. In the latter we have found that the Hawking radiation affects, although very slightly, the oscillations for neutrinos propagating from the asymptotic past infinity to the asymptotic future infinity. As a general result, it is found that when all the quantum field theoretical effects on neutrino mixing can be neglected, the gravitational background only affects the phase of the oscillations, constistently with previous analyses.