# The flavor vacuum in the expanding universe and dark matter

Antonio Capolupo<EMAIL_ADDRESS>Dipartimento di Fisica “E.R. Caianiello”
Università di Salerno, and INFN – Gruppo Collegato di Salerno, Via Giovanni
Paolo II, 132, 84084 Fisciano (SA), Italy  Sante Carloni
<EMAIL_ADDRESS>DIME Sezione Metodi e Modelli Matematici, Università di
Genova, Via All’Opera Pia 15, 16145 Genova, Italy  Aniello Quaranta
<EMAIL_ADDRESS>Dipartimento di Fisica “E.R. Caianiello” Università di
Salerno, and INFN – Gruppo Collegato di Salerno, Via Giovanni Paolo II, 132,
84084 Fisciano (SA), Italy

###### Abstract

We analyze fermion mixing in the framework of field quantization in curved
spacetime. We compute the expectation value of the energy momentum tensor of
mixed fermions on the flavor vacuum. We consider spatially flat Friedmann-
Lemaitre-Robertson-Walker metrics, and we show that the energy–momentum tensor
of the flavor vacuum is diagonal and conserved. Therefore it can be
interpreted as the effective energy–momentum tensor of a perfect fluid. In
particular, assuming a fixed De Sitter background, the equation of state of
the fluid is consistent with that of dust and cold dark matter. Our results
establish a new link between quantum effects and classical fluids, and
indicate that the flavor vacuum of mixed fermions may represent a new
component of dark matter.

##  I Introduction

Among the open issues of modern cosmology is the understanding of the dark
components of the universe. Most of the total energy density is shared between
“dark energy”, that drives the accelerated expansion of the universe [1]-[7],
and “dark matter”, the non–baryonic matter which is responsible for holding
galaxies and clusters together [8]-[25]. A conclusive explanation for this
“dark universe” is still lacking. Several proposals, often very different in
nature and perspective, have been put forward to account for dark matter. New
particles such as axions, axionlike particles [26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36] and supersymmetric partners [37], arising from extensions of the
standard model, might explain the missing matter. Another possible dark matter
component is represented by massive compact objects like primordial black
holes [38]. Other proposals rely on pure quantum field theoretical effects and
the non–trivial structure of the vacuum of flavor fermion mixing [39, 40, 41,
42, 43]. The relevance of fermion fields in astrophysical and cosmological
contexts is well represented by the role of neutrinos. These particles are a
valuable source of information as they are expected to play an important role
in investigating astrophysical processes. They can be used as a test of the
standard cosmological model [45, 46] even at early time [44, 45] and have been
recognized as a possible component of dark energy [47]. Moreover the
corresponding vacuum has been identified as a possible dark matter component
in flat space [39, 40]. Neutrinos are known to oscillate among their three
flavors and their general oscillation formulas in curved space-time have been
derived in refs. [48, 52, 49, 50, 51] together with a number of other aspects
of this phenomenon.

In the following, starting from the analysis of the possible contribution of
the flavor vacuum to the dark matter, previously studied in Minkowski
spacetime [39, 40], we study the behaviour of the flavor vacuum in the case of
curved background. We compute, in the context of an homogeneous and isotropic
spacetime, represented by a spatially flat Friedmann–Lemaitre–Robertson–Walker
metric, the expectation value of the energy momentum tensor of neutrinos on
the flavor (mixed) vacuum. We show that this expectation value is a diagonal
tensor and satisfies the Bianchi identities. Consequently it behaves as an
effective stress energy tensor akin to that of a perfect fluid and can enter
the Einstein equations as a regular source term. In particular, assuming a
fixed De Sitter background, the equation of state of such a fluid is that of
dust or cold dark matter ($w=0$).

These results are not a mere generalization of the studies conducted in flat
space [39, 40, 41, 42, 43], since the underlying quantum field theory of
fermion mixing is much more involved, with respect to flat space, and the
properties of the flavor vacuum depend critically on the curved background
considered. The computation performed on the De Sitter background already
gives an indication that the flavor vacuum may play a role within the dark
sector of the universe. Indeed the flavor vacuum brings along an additional
energy density, and, being pressure–less, may be associated with a dark matter
component. According to our results, a possible constituent of dark matter is
then represented by the energy of the vacuum of mixed fields. In future works
we will analyze the energy–momentum tensor associated to the fermion flavor
vacuum in spacetimes that describe large scale structures, like spiral
galaxies. In that context the fluid related to the flavor vacuum can be more
properly identified with a dark matter component.

The paper is organized as following. Section II contains a summary of the
properties of the quantum Dirac equation in curved spacetime. Section III
contains the field quantization of the flavor fields and the introduction of
the oscillation components. In Section IV the component of the expectation
value of the quantum stress energy tensor is derived. Section V contains the
application to a specific background metric i.e. the de Sitter spacetime and
the derivation of the corresponding exact solution for the component of the
stress energy tensor, together with some consideration on its regularization.
Finally Section VI is dedicated to the conclusions. In the following, we will
use the $+---$ signature, and we will assume that greek indices run from
$0,...,3$ lower case Latin indices from $1,...,3$ and upper case Latin indices
from $1,...,4$. These last indices will represent tetrad indices. The symbol
$[,]$ denotes the matrix commutator so that, e.g.,
$[\gamma_{A},\gamma_{B}]=\gamma_{A}\gamma_{B}-\gamma_{B}\gamma_{A}$.

##  II The Dirac Equation and its solutions in flat FLRW spacetime

For the reader’s convenience we start this section by setting the notation and
introducing the metric of interest. We then elaborate on the Dirac equation
and discuss the general form of the solution. We will focus on the spatially
flat Friedmann-Lemaître-Robertson-Walker spacetime (we use the $+---$
signature) described by the metric
$g_{\mu\nu}=\mathrm{diag}\left(1,-C^{2}(t),-C^{2}(t),-C^{2}(t)\right).$ where
$C(t)$ is the scale factor. For our purposes it will be useful to express this
metric in terms of the conformal time $\tau$ defined as
$d\tau=\frac{dt}{C(t)}$ with range $\tau\in(-\infty,\infty)$ corresponding to
$t\in(-\infty,\infty)$. Using $\tau$, the line element reads

$ds^{2}=C^{2}(\tau)[d\tau^{2}-dx^{2}-dy^{2}-dz^{2}]$ (1)

when expressed in cartesian spatial coordinates and in conformal time $\tau$.

In order to write the Dirac equation, we need to choose a tetrad field
$e^{A}_{\mu}$ with the property $g_{\mu\nu}=e^{A}_{\mu}e^{B}_{\nu}\eta_{AB}$
Given the metric of Eq.(1) a convenient choice of tetrads is

$e^{A}_{\mu}=C(\tau)\delta^{A}_{\mu}\ ,$ (2)

where the Kronecker symbol $\delta^{A}_{\mu}$ signals that the only non-zero
component of $e^{A}$ is the one for which the Lorentz index $A$ and the
spacetime index $\mu$ coincide.

Using the tetrads above one can define the generalization of gamma matrices
$\gamma^{A}$ to the case of curved spacetime
$\tilde{\gamma}^{\mu}=e^{\mu}_{A}\gamma^{A}$

Finally, in order to write the Dirac equation, we also need the spin
connections $\omega_{\mu}^{AB}=e^{A}_{\nu}\Gamma^{\nu}_{\sigma\mu}e^{\sigma
B}+e^{A}_{\nu}\partial_{\mu}e^{\nu B}\ .$ The spinorial covariant derivative
is defined with the aid of the spin connections as
$D_{\mu}\psi=\partial_{\mu}\psi+\Gamma_{\mu}\psi\,,\ \ \ \ \ \ \
\Gamma_{\mu}=\frac{1}{8}\omega_{\mu}^{AB}[\gamma_{A},\gamma_{B}].$ We can now
write down the Dirac equation

$i\tilde{\gamma}^{\mu}(x)D_{\mu}\psi-m\psi=0\ .$ (3)

The above equation can be generated by the Lagrangian

$\mathcal{L}=\sqrt{-g}\left\\{\frac{i}{2}\left[\bar{\psi}\tilde{\gamma}^{\mu}(x)D_{\mu}\psi-
D_{\mu}\bar{\psi}\tilde{\gamma}^{\mu}(x)\psi\right]-m\bar{\psi}\psi\right\\}$
(4)

and we can define the energy-momentum tensor of the spinor field as the
variation of the above Lagrangian with respect to $\psi$ and $\bar{\psi}$

$T_{\mu\nu}=\frac{i}{2}\left\\{\bar{\psi}\tilde{\gamma}_{\mu}(x)D_{\nu}\psi+\bar{\psi}\tilde{\gamma}_{\nu}(x)D_{\mu}\psi-
D_{\mu}\bar{\psi}\tilde{\gamma}_{\nu}(x)\psi-
D_{\nu}\bar{\psi}\tilde{\gamma}_{\mu}(x)\psi\right\\}\ .$ (5)

In the metric (1) the Dirac equation reads

$\left(i\gamma^{0}\partial_{\tau}+\frac{3i}{2}\frac{\partial_{\tau}C}{C}\gamma^{0}+i\gamma^{j}\partial_{j}-mC\right)\psi=0\
.$ (6)

We remark that the Dirac equation (6) holds for the quantum field theoretic
free Dirac field of mass $m$. At odds with the _classical_ Dirac field, we
cannot drop the spatial derivatives by arguing that space-dependent quantities
cannot enter the right hand side of the Einstein equation if the metric
depends only on time. Indeed _all_ the solutions to Eq. (6) must be considered
for the field expansion, including those which depend explicitly on the
spatial coordinates. Consistency with the (time-dependent only) metric is then
achieved for the expectation value of the physical observables, including the
energy-momentum tensor, which turn out to depend only upon time.

The spatial dependence of this equation suggests that we look for plane wave
solutions $\psi\propto e^{i\boldsymbol{p}\cdot\boldsymbol{x}}$, where the mode
label $\boldsymbol{p}$ is a $3$-vector that can be thought as the would-be
plane wave momentum when $C(\tau)$ reduces to a constant, and
”$\boldsymbol{a}\cdot\boldsymbol{b}$” is a shorthand notation for
$\sum_{j=1,2,3}a_{j}b_{j}$. We remark that the actual momentum that is
istanteneously carried by the mode with label $\boldsymbol{p}$ is the comoving
momentum111This is most conveniently seen by inserting the plane wave ansatz
in Eq. (6) and reverting to the coordinate time $t$:
$\left[i\gamma^{0}\left(\partial_{t}+\frac{3i}{2}\frac{\partial_{t}C}{C}\right)-\frac{\gamma^{j}p_{j}}{C}-m\right]\psi_{\boldsymbol{p}}=0\
.$ In this form the equation resembles the Dirac equation in flat space, with
a time-dependent potential $\frac{3i}{2}\frac{\partial_{t}C}{C}$ and with the
quantity $\frac{p_{j}}{C}$ playing the role of an instantaneous momentum.
$\frac{\boldsymbol{p}}{C(\tau)}$.

Using the helicity eigenbispinors $\xi_{\lambda}$ defined as

$\frac{\boldsymbol{\sigma}\cdot\boldsymbol{p}}{p}\xi_{\lambda}=\lambda\xi_{\lambda}\
,$ (7)

the solution of (6) can be written as the combination of the positive and
negative energy solutions in the form

$u_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{x})=e^{i\boldsymbol{p}\cdot\boldsymbol{x}}\begin{pmatrix}f_{p}(\tau)\xi_{\lambda}(\hat{p})\\\
g_{p}(\tau)\lambda\xi_{\lambda}(\hat{p})\end{pmatrix}\ \ \ \
v_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{x})=e^{i\boldsymbol{p}\cdot\boldsymbol{x}}\begin{pmatrix}g^{*}_{p}(\tau)\xi_{\lambda}(\hat{p})\\\
-f^{*}_{p}(\tau)\lambda\xi_{\lambda}(\hat{p})\end{pmatrix}\ .$ (8)

Here $\hat{p}=\frac{\boldsymbol{p}}{p}$ denotes the unit vector in the
direction of $\boldsymbol{p}$, $\boldsymbol{\sigma}$ is the vector of Pauli
matrices and $\lambda=\pm 1$. Notice that for each $\lambda=\pm 1$, the
quantity $\xi^{\dagger}_{\lambda}(\hat{p})\sigma_{i}\xi_{\lambda}(\hat{p})$ is
an odd function of $p_{i}$, for $i=1,2,3$. In particular,
$\xi^{\dagger}_{\lambda}(\hat{p})\sigma_{i}\xi_{\lambda}(\hat{p})$ changes
sign when the momentum is reversed $\boldsymbol{p}\rightarrow-\boldsymbol{p}$.
The statement can be proven by direct calculation, as shown in Appendix B.
Inserting Eq.(8) in Eq.(6), and using the defining property of the helicity
bispinors, we can write (6) as the system

$\begin{split}i\partial_{\tau}f_{\boldsymbol{p}}&=\left(mC-\frac{3i}{2}\frac{\partial_{\tau}C}{C}\right)f_{\boldsymbol{p}}+pg_{\boldsymbol{p}}\\\
i\partial_{\tau}g_{\boldsymbol{p}}&=\left(-mC-\frac{3i}{2}\frac{\partial_{\tau}C}{C}\right)g_{\boldsymbol{p}}+pf_{\boldsymbol{p}}\
,\\\
i\partial_{\tau}f^{*}_{p}&=\left(-mC-\frac{3i}{2}\frac{\partial_{\tau}C}{C}\right)f^{*}_{p}-pg^{*}_{p}\\\
i\partial_{\tau}g^{*}_{p}&=\left(mC-\frac{3i}{2}\frac{\partial_{\tau}C}{C}\right)g^{*}_{p}-pf^{*}_{p}\
.\end{split}$ (9)

The above equations show that the functions $f,g$ depend only on the modulus
$p$ of the $3$-momentum. Thus, from now on, we will drop the vector index
$\boldsymbol{p}$ in favor of the scalar index $p$.

The scalar product between two solutions of the Dirac equation $A,B$ is
defined as

$(A,B)_{\tau}=\int_{\Sigma_{\tau}}d^{3}x\sqrt{-g}\bar{A}\tilde{\gamma}^{\tau}(x)B\
,$ (10)

where the integration is to be carried over a hypersurface $\Sigma_{\tau}$ of
constant conformal time $\tau$. If $A$ and $B$ are solutions of the same Dirac
equation, the scalar product $(A,B)_{\tau}$ is independent of $\tau$. This is
no longer true if $A$ and $B$ are solutions of distinct Dirac equations.
Taking into account that the determinant is $g=-C^{8}$ and adopting the
normalization

$|f_{p}|^{2}+|g_{p}|^{2}=\frac{1}{\left(2\pi C\right)^{3}}$ (11)

we obtain the following orthonormality and completeness relations,
respectively (for details, see Appendix A):

$\displaystyle(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\delta_{\lambda,\lambda^{\prime}}\
,\qquad(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=(v_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=0\
,\qquad(v_{\boldsymbol{p},\lambda},v_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\delta_{\lambda,\lambda^{\prime}}\
,$ (12)
$\displaystyle\sum_{\lambda}\left(u_{\boldsymbol{p},\lambda}u^{\dagger}_{\boldsymbol{p},\lambda}+v_{\boldsymbol{p},\lambda}v^{\dagger}_{\boldsymbol{p},\lambda}\right)$
$\displaystyle=$ $\displaystyle\frac{1}{\left(2\pi
C\right)^{3}}\begin{pmatrix}\mathbb{I}&0\\\ 0&\mathbb{I}\end{pmatrix}\ .$ (13)

Eq.(13) shows that the set of solutions (8) provides a basis not only for the
solution space of the Dirac equation, but also for the space of $4$-spinors.

The system of equations (9) can be further simplified by introducing the
functions

$\phi_{p}=C^{\frac{3}{2}}f_{p}\ \ \ \ \gamma_{p}=C^{\frac{3}{2}}g_{p}$ (14)

since then the system becomes

$\begin{split}\partial_{\tau}\phi_{p}&=-imC\phi_{p}-ip\gamma_{p}\\\
\partial_{\tau}\gamma_{p}&=imC\gamma_{p}-ip\phi_{p}\\\
\partial_{\tau}\phi^{*}_{p}&=imC\phi^{*}_{p}+ip\gamma^{*}_{p}\\\
\partial_{\tau}\gamma^{*}_{p}&=-imC\gamma^{*}_{p}+ip\phi^{*}_{p}\end{split}$
(15)

and the normalization condition (11) becomes

$|\phi_{p}|^{2}+|\gamma_{p}|^{2}=\frac{1}{(2\pi)^{3}}\ .$ (16)

The first two of Eqs.(15) can be combined to give two second order equation
for $\phi_{p}$

$\partial_{\tau}^{2}\phi_{p}+\left(im\partial_{\tau}C+p^{2}+m^{2}C^{2}\right)\phi_{p}=0\
.$ (17)

In the same way, we can obtain a second-order equation for $\phi^{*}_{p}$.

We conclude the section by introducing a bilinear form of the solutions of the
Dirac Equation, which is convenient for the study of the energy-momentum
tensor. Given two solutions $A,B$ of the Dirac equation, we define the
auxiliary tensor as

${}L_{\mu\nu}(A,B)=\bar{A}\tilde{\gamma}_{\mu}(x)D_{\nu}B+\bar{A}\tilde{\gamma}_{\nu}(x)D_{\mu}B-D_{\mu}\bar{A}\tilde{\gamma}_{\nu}(x)B-D_{\nu}\bar{A}\tilde{\gamma}_{\mu}(x)B\
.$ (18)

The properties of the auxiliary tensor are analyzed in Appendix C.

##  III Quantization of flavor fields

In the following, we start by quantizing a single Dirac field of definite
mass, and then we report the analysis of the main properties of two flavor
mixed fields. For more details, we refer to [48].

###  III.1 Dirac field

Since the set of solutions (8)
$\\{u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda}\\}_{\boldsymbol{p},\lambda}$
is complete, any solution of the (linear) Dirac equation can be written as a
linear combination of these modes. In particular, this is true for the Dirac
field: $\psi(x)=\sum_{\lambda}\int
d^{3}p\left(A_{\boldsymbol{p,\lambda}}u_{\boldsymbol{p},\lambda}+B^{*}_{-\boldsymbol{p},\lambda}v_{\boldsymbol{p},\lambda}\right)$.
In this equation, the notation $B^{*}_{-\boldsymbol{p},\lambda}$ is chosen to
remark that $v_{\boldsymbol{p},\lambda}$ describes an antiparticle with
momentum $-\boldsymbol{p}$. Notice that all the spacetime dependence is in the
modes, while the coefficients are independent of both space and time
coordinates.

Quantization is achieved, as usual, by promoting the field, and thus the
expansion coefficients to operators

$\psi(x)=\sum_{\lambda}\int
d^{3}p\left(A_{\boldsymbol{p,\lambda}}u_{\boldsymbol{p},\lambda}+B^{\dagger}_{-\boldsymbol{p},\lambda}v_{\boldsymbol{p},\lambda}\right)\
.$ (19)

and imposing the canonical anticommutation relations. The momentum conjugate
to $\psi(x)$, according to the Lagrangian (4), is
$\pi_{\psi}(x)=iC^{3}\psi^{\dagger}(x)\ ,$ so that the canonical
anticommutation relations to be imposed are

$\\{\psi_{A}(\tau,\boldsymbol{x}),\pi_{\psi
B}(\tau,\boldsymbol{x^{\prime}})\\}=iC^{3}\\{\psi_{A}(\tau,\boldsymbol{x}),\psi^{\dagger}_{B}(\tau,\boldsymbol{x^{\prime}})\\}=i\delta_{AB}\delta^{3}(\boldsymbol{x}-\boldsymbol{x^{\prime}})\
.$ (20)

Here the indices $A,B$ are referred to the spinor components $A,B=1,2,3,4$. It
is easy to show that the relations (20) are satisfied if one imposes the
following anticommutation relations on the coefficients:

$\\{A_{\boldsymbol{p,\lambda}},A^{\dagger}_{\boldsymbol{q},\lambda^{\prime}}\\}=\\{B_{\boldsymbol{p,\lambda}},B^{\dagger}_{\boldsymbol{q},\lambda^{\prime}}\\}=\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})$
(21)

with all the other anticommutators vanishing (see Appendix …). The field
expansion (19) defines the vacuum state $\ket{0}$ as the state annihilated by
all the annihilation operators
$A_{\boldsymbol{p},\lambda}\ket{0}=B_{\boldsymbol{p},\lambda}\ket{0}=0\,,\ \ \
\ \forall\boldsymbol{p},\lambda\ .$ It is important to stress that the
definition of the vacuum state depends critically on the choice of the field
modes. Specifically, its particle interpretation is tied to the boundary
conditions specified on the solutions (8). Another kind of field expansion is
possible if one assumes a specific time evolution of the modes, as is done
within the adiabatic approximation (see, e.g., [53, 54, 55, 56]). Contrary to
our field expansion (19), the annihilation operators (and thus the vacuum) are
thereby endowed with a specific time dependence. The expansion (19), instead,
does not assume any particular time dependence and is, therefore, more
general. The two expansions can be made to coincide at a given time. The
quantized energy-momentum tensor is obtained by inserting the field expansion
in the definition (5):

$\displaystyle T_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{i}{2}\sum_{\lambda,\lambda^{\prime}}\int d^{3}p\int
d^{3}q\\{A^{\dagger}_{\boldsymbol{p},\lambda}A_{\boldsymbol{q},\lambda^{\prime}}L_{\mu\nu}(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})+A^{\dagger}_{\boldsymbol{p},\lambda}B^{\dagger}_{-\boldsymbol{q},\lambda^{\prime}}L_{\mu\nu}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{q},\lambda^{\prime}})$
(22)
$\displaystyle+B_{-\boldsymbol{p},\lambda}A_{\boldsymbol{q},\lambda^{\prime}}L_{\mu\nu}(v_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})+B_{-\boldsymbol{p},\lambda}B^{\dagger}_{-\boldsymbol{q},\lambda^{\prime}}L_{\mu\nu}(v_{\boldsymbol{p},\lambda},v_{\boldsymbol{q},\lambda^{\prime}})\\}\
,$

where we have used the definition (18) of $L_{\mu\nu}$. The Hamiltonian
density, defined with respect to $\partial_{\tau}$ corresponds to the
component $T_{\tau}^{\tau}$ of the above equation. It follows immediately that
the vacuum is not an eigenstate of the Hamiltonian unless
$L^{\tau}_{\tau}(v_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})=0$,
due to the
$A^{\dagger}_{\boldsymbol{p},\lambda}B^{\dagger}_{-\boldsymbol{q},\lambda^{\prime}}$
term. In particular, the vacuum is unstable under the creation of particle-
antiparticle pairs with opposite momentum, a result known from the previous
analyses. This is a reflection of the non-invariance under time translations.
On the other hand, from the relation: $L_{\tau
i}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})=0$, derived in
Appendix C (cfr. Eq.(94)), it is clear that the vacuum respects the residual
translational symmetry in the spatial coordinates, i.e.
$\boldsymbol{P}\ket{0}=0$.

###  III.2 The flavor fields

Up to now, our considerations have been restricted to a single Dirac field of
definite mass. To introduce the flavor fields, we follow [48] and start by
introducing two Dirac fields with distinct masses $m_{1},m_{2}$. We consider
two flavors for simplicity but the analysis can be easily generalized to three
flavors. The theory of two free massive Dirac fields is just the product of
two copies of the theory for a single Dirac field. All the relations discussed
above remain valid, provided we assign a new index $j=1,2$ to all the
quantities involved
$\psi_{j},m_{j},A_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j}$ and
so on. All the previous relations are index-wise valid for $j=1,2$. The mass
vacuum that we now denote explicitly as $\ket{0_{M}}$ is annihilated by all
the annihilators for each index $j$. The total energy-momentum tensor is
simply the sum of two copies of Eq.(22), one for each $j$. We also require
that each mode of the field $1$ is related to the corresponding mode of the
field $2$, with the same labels $\boldsymbol{p},\lambda$ by the substitution
$m_{1}\rightarrow m_{2}$, and vice-versa.

The flavor fields are then introduced via the rotation

$\displaystyle\psi_{e}(x)$ $\displaystyle=$
$\displaystyle\cos(\theta)\psi_{1}(x)+\sin(\theta)\psi_{2}(x)$ $\displaystyle\
\psi_{\mu}(x)$ $\displaystyle=$
$\displaystyle\cos(\theta)\psi_{2}(x)-\sin(\theta)\psi_{1}(x)\ ,$ (23)

where $\theta$ is the two-flavor mixing angle. At the quantum level, the
rotation is employed by the mixing generator

$\mathcal{I}_{\theta}(\tau)=exp\left\\{\theta\left[(\psi_{1},\psi_{2})_{\tau}-(\psi_{2},\psi_{1})_{\tau}\right]\right\\}\
.$ (24)

Here $(\psi_{2},\psi_{1})_{\tau}$ stands for the scalar product at the $\tau$
hypersurface (recall that for fields of distinct masses the product _does_
depend on $\tau$). The flavor fields are then

$\displaystyle\psi_{e}(x)$ $\displaystyle=$
$\displaystyle\mathcal{I}^{-1}_{\theta}\psi_{1}(x)\mathcal{I}_{\theta}$
$\displaystyle\psi_{\mu}(x)$ $\displaystyle=$
$\displaystyle\mathcal{I}^{-1}_{\theta}\psi_{2}(x)\mathcal{I}_{\theta}\ .$
(25)

The action of the generator also defines the flavor annihilators
($A_{\boldsymbol{p},\lambda;e}=\mathcal{I}^{-1}_{\theta}A_{\boldsymbol{p},\lambda;1}\mathcal{I}_{\theta}$
and similar) and the flavor vacuum, as

$\ket{0_{F}}=\mathcal{I}^{-1}_{\theta}\ket{0_{M}}\ .$ (26)

Notice that, contrary to the mass vacuum, the flavor vacuum has an explicit
$\tau$ dependence. The terminology ’flavor vacuum’ is justified in that this
state is annihilated by all the flavor annihilators.

##  IV VEV of the energy-momentum tensor on the flavor vacuum

We are specifically interested in the contributions that flavor mixing induces
on the energy-momentum tensor. More precisely, we ask what is the expectation
value of the energy-momentum tensor on the state corresponding to the absence
of flavor neutrinos at a given time $\ket{O_{F}(\tau_{0})}$. We stress that in
this calculation, we assume a fixed but arbitrary expansion of the mass
fields, and therefore a fixed but arbitrary choice of the mass vacuum. The
effect of a change in the mass representation on the flavor fields is known
[48], and once the result is computed for a given representation, one can
implement the adequate transformations to get the result in other mass
representations. Likewise, we keep the time $\tau_{0}$ arbitrary and distinct
from the time argument $\tau$ of the energy-momentum tensor. The quantity we
wish to compute is then

$\mathbb{T}_{\mu\nu}=\bra{0_{F}(\tau_{0})}T_{\mu\nu}\ket{0_{F}(\tau_{0})}\ ,$
(27)

where $T_{\mu\nu}$ is given by Eq.(22). We remark that the only sensible
definition for the energy-momentum tensor is the one (22) in terms of the
fields with definite mass. Let us analyze the typical term in equation (27).
It has the form

$\bra{0_{F}(\tau_{0})}A^{\dagger}_{\boldsymbol{p},\lambda;1}A_{\boldsymbol{q},\lambda^{\prime};1}\ket{0_{F}(\tau_{0})}\
.$ (28)

Using the definition of the flavor vacuum this equals

$\displaystyle\bra{0_{M}}\mathcal{I}_{\theta}(\tau_{0})A^{\dagger}_{\boldsymbol{p},\lambda;1}A_{\boldsymbol{q},\lambda^{\prime};1}\mathcal{I}^{-1}_{\theta}(\tau_{0})\ket{0_{M}}$
$\displaystyle=$
$\displaystyle\bra{0_{M}}\mathcal{I}_{\theta}(\tau_{0})A^{\dagger}_{\boldsymbol{p},\lambda;1}\mathcal{I}^{-1}_{\theta}(\tau_{0})\mathcal{I}_{\theta}(\tau_{0})A_{\boldsymbol{q},\lambda^{\prime};1}\mathcal{I}^{-1}_{\theta}(\tau_{0})\ket{0_{M}}$
(29) $\displaystyle=$
$\displaystyle\bra{0_{M}}\mathcal{I}^{-1}_{-\theta}(\tau_{0})A^{\dagger}_{\boldsymbol{p},\lambda;1}\mathcal{I}_{-\theta}(\tau_{0})\mathcal{I}^{-1}_{-\theta}(\tau_{0})A_{\boldsymbol{q},\lambda^{\prime};1}\mathcal{I}_{-\theta}(\tau_{0})\ket{0_{M}}$

where we have used that (eq. 24)
$\mathcal{I}^{-1}_{\theta}=\mathcal{I}_{-\theta}$. Now the operator
$\mathcal{I}^{-1}_{-\theta}(\tau_{0})A_{\boldsymbol{p},\lambda;1}\mathcal{I}_{-\theta}(\tau_{0})$
is just the mass annihilator transformed according to a mixing transformation
with angle $-\theta$. Knowing the transformation rule in terms of $\theta$
[48] we can easily write down the transformed operators:

$\begin{split}\mathcal{I}_{-\theta}^{-1}(\tau_{0})A_{\boldsymbol{p},\lambda;1}\mathcal{I}_{-\theta}(\tau_{0})&=\cos(\theta)A_{\boldsymbol{p},\lambda;1}-\sin(\theta)\left(\Lambda^{*}_{\boldsymbol{p}}(\tau_{0})A_{\boldsymbol{p},\lambda;2}+\Xi_{\boldsymbol{p}}(\tau_{0})B^{\dagger}_{-\boldsymbol{p},\lambda;2}\right)\\\
\mathcal{I}_{-\theta}^{-1}(\tau_{0})A_{\boldsymbol{p},\lambda;2}\mathcal{I}_{-\theta}(\tau_{0})&=\cos(\theta)A_{\boldsymbol{p},\lambda;2}+\sin(\theta)\left(\Lambda_{\boldsymbol{p}}(\tau_{0})A_{\boldsymbol{p},\lambda;1}-\Xi_{\boldsymbol{p}}(\tau_{0})B^{\dagger}_{-\boldsymbol{p},\lambda;1}\right)\\\
\mathcal{I}_{-\theta}^{-1}(\tau_{0})B_{-\boldsymbol{p},\lambda;1}\mathcal{I}_{-\theta}(\tau_{0})&=\cos(\theta)B_{-\boldsymbol{p},\lambda;1}-\sin(\theta)\left(\Lambda^{*}_{\boldsymbol{p}}(\tau_{0})B_{-\boldsymbol{p},\lambda;2}-\Xi_{\boldsymbol{p}}(\tau_{0})A^{\dagger}_{\boldsymbol{p},\lambda;2}\right)\\\
\mathcal{I}_{-\theta}^{-1}(\tau_{0})B_{-\boldsymbol{p},\lambda;2}\mathcal{I}_{-\theta}(\tau_{0})&=\cos(\theta)B_{-\boldsymbol{p},\lambda;2}+\sin(\theta)\left(\Lambda_{\boldsymbol{p}}(\tau_{0})B_{-\boldsymbol{p},\lambda;1}+\Xi_{\boldsymbol{p}}(\tau_{0})A^{\dagger}_{\boldsymbol{p},\lambda;1}\right)\end{split}$
(30)

while the transformation rule for the adjoint operators can be obtained by
considering the adjoint equations. The Bogoliubov coefficients are defined as
the inner products

$\begin{split}\delta^{3}(\boldsymbol{0})\Lambda_{\boldsymbol{p}}(\tau)&=\left(u_{\boldsymbol{p},\lambda;2},u_{\boldsymbol{p},\lambda;1}\right)_{\tau}=\left(v_{\boldsymbol{p},\lambda;1},v_{\boldsymbol{p},\lambda;2}\right)_{\tau}\\\
\delta^{3}(\boldsymbol{0})\Xi_{\boldsymbol{p}}(\tau)&=\left(u_{\boldsymbol{p},\lambda;1},v_{\boldsymbol{p},\lambda;2}\right)_{\tau}=-\left(u_{\boldsymbol{p},\lambda;2},v_{\boldsymbol{p},\lambda;1}\right)_{\tau}\
,\end{split}$ (31)

where the notation anticipates that they don’t depend on $\lambda$ (we will
see that they actually depend only on the modulus of the momentum $p$). The
$\delta^{3}(\boldsymbol{0})$ factor is a reminiscence of the more general
expressions involving distinct momenta
$\Lambda_{\boldsymbol{p},\boldsymbol{q}}\propto\delta^{3}(\boldsymbol{p}-\boldsymbol{q})$
[48]. The delta factor is absorbed by a corresponding momentum integration in
the Eqs.(30), leaving only the finite coefficients defined via Eq.(31). These
coefficients satisfy
$|\Lambda_{\boldsymbol{p}}|^{2}+|\Xi_{\boldsymbol{p}}|^{2}=1$ for all
$\boldsymbol{p},\tau$. With the aid of Eqs.(30), we can evaluate all the
expectation values appearing in Eq.(27):

$\begin{split}\bra{0_{F}(\tau_{0})}A^{\dagger}_{\boldsymbol{p},\lambda;j}A_{\boldsymbol{q};\lambda^{\prime};j}\ket{0_{F}(\tau_{0})}&=\sin^{2}\theta|\Xi_{\boldsymbol{p}}(\tau_{0})|^{2}\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\
\,,\ \ \ \forall j\\\
\bra{0_{F}(\tau_{0})}B^{\dagger}_{-\boldsymbol{p},\lambda;j}B_{-\boldsymbol{q};\lambda^{\prime};j}\ket{0_{F}(\tau_{0})}&=\sin^{2}\theta|\Xi_{\boldsymbol{p}}(\tau_{0})|^{2}\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\
\,,\ \ \ \forall j\\\
\bra{0_{F}(\tau_{0})}A^{\dagger}_{\boldsymbol{p},\lambda;1}B^{\dagger}_{-\boldsymbol{q};\lambda^{\prime};1}\ket{0_{F}(\tau_{0})}&=\sin^{2}\theta\Xi^{*}_{\boldsymbol{p}}(\tau_{0})\Lambda_{\boldsymbol{p}}(\tau_{0})\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\\\
\bra{0_{F}(\tau_{0})}A^{\dagger}_{\boldsymbol{p},\lambda;2}B^{\dagger}_{-\boldsymbol{q};\lambda^{\prime};2}\ket{0_{F}(\tau_{0})}&=-\sin^{2}\theta\Xi^{*}_{\boldsymbol{p}}(\tau_{0})\Lambda^{*}_{\boldsymbol{p}}(\tau_{0})\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\\\
\bra{0_{F}(\tau_{0})}B_{-\boldsymbol{p},\lambda;1}A_{\boldsymbol{q};\lambda^{\prime};1}\ket{0_{F}(\tau_{0})}&=\sin^{2}\theta\Xi_{\boldsymbol{p}}(\tau_{0})\Lambda^{*}_{\boldsymbol{p}}(\tau_{0})\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\\\
\bra{0_{F}(\tau_{0})}B_{-\boldsymbol{p},\lambda;2}A_{\boldsymbol{q};\lambda^{\prime};2}\ket{0_{F}(\tau_{0})}&=-\sin^{2}\theta\Xi_{\boldsymbol{p}}(\tau_{0})\Lambda_{\boldsymbol{p}}(\tau_{0})\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\
.\\\ \end{split}$ (32)

It is convenient to give the result by splitting out the pure mixing
component:

$\displaystyle\mathbb{T}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\mathbb{T}_{\mu\nu}^{(MIX)}+\mathbb{T}_{\mu\nu}^{(N)}$ (33)
$\displaystyle\mathbb{T}_{\mu\nu}^{(MIX)}$ $\displaystyle=$
$\displaystyle\frac{i}{2}\sin^{2}\theta\,\sum_{\lambda}\int
d^{3}p\\{|\Xi_{\boldsymbol{p}}(\tau_{0})|^{2}\sum_{j=1,2}\left(L_{\mu\nu}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})-L_{\mu\nu}(v_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})\right)$
(34) $\displaystyle+$ $\displaystyle\
\Xi^{*}_{\boldsymbol{p}}(\tau_{0})\Lambda_{\boldsymbol{p}}(\tau_{0})L_{\mu\nu}(u_{\boldsymbol{p},\lambda;1},v_{\boldsymbol{p},\lambda;1})+\Xi_{\boldsymbol{p}}(\tau_{0})\Lambda^{*}_{\boldsymbol{p}}(\tau_{0})L_{\mu\nu}(v_{\boldsymbol{p},\lambda;1},u_{\boldsymbol{p},\lambda;1})$
$\displaystyle-$
$\displaystyle\Xi^{*}_{\boldsymbol{p}}(\tau_{0})\Lambda^{*}_{\boldsymbol{p}}(\tau_{0})L_{\mu\nu}(u_{\boldsymbol{p},\lambda;2},v_{\boldsymbol{p},\lambda;2})-\Xi_{\boldsymbol{p}}(\tau_{0})\Lambda_{\boldsymbol{p}}(\tau_{0})L_{\mu\nu}(v_{\boldsymbol{p},\lambda;2},u_{\boldsymbol{p},\lambda;2})\\}$
$\displaystyle\mathbb{T}_{\mu\nu}^{(N)}$ $\displaystyle=$
$\displaystyle\frac{i}{2}\sum_{\lambda}\sum_{j=1,2}\int
d^{3}pL_{\mu\nu}(v_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})\
.$ (35)

The last term comes from applying the anticommutator relation
$\left\\{B^{\dagger}_{-\boldsymbol{p},\lambda;j},B_{-\boldsymbol{q},\lambda^{\prime};j}\right\\}=\delta_{\lambda\lambda^{\prime}}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})$
to the $BB^{\dagger}$ term. While the remaining terms, regrouped under the
symbol $\mathbb{T}_{\mu\nu}^{(MIX)}$ show an explicit dependence on
$\sin^{2}\theta$, and are therefore zero in absence of mixing, the last term
is present independently of mixing. Indeed, it is easy to check that
$\mathbb{T}^{(N)}_{\mu\nu}$ represents the expectation value of the energy-
momentum tensor on the _mass_ vacuum:

$\mathbb{T}^{(N)}_{\mu\nu}=\bra{0_{M}}T_{\mu\nu}\ket{0_{M}}\ .$ (36)

The origin of this term is an ordering ambiguity in the energy-momentum tensor
of the quantized fields. To understand its significance, let us consider the
flat space limit, where $v_{\boldsymbol{p},\lambda;j}\propto
e^{i\omega_{p;j}t}$ with $\omega_{p;j}=\sqrt{p^{2}+m_{j}^{2}}$. The auxiliary
tensor becomes

$L_{\mu\nu}(A,B)=\bar{A}\gamma_{\mu}\partial_{\nu}B+\bar{A}\gamma_{\nu}\partial_{\mu}B-\partial_{\mu}\bar{A}\gamma_{\nu}B-\partial_{\nu}\bar{A}\gamma_{\mu}B$
(37)

where both the Dirac matrices and the derivatives are the ordinary flat space
ones. In particular one has

$L_{00}(v_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})=2(v_{\boldsymbol{p},\lambda;j}^{\dagger}\partial_{t}v_{\boldsymbol{p},\lambda;j}-\partial_{t}v_{\boldsymbol{p},\lambda;j}^{\dagger}v_{\boldsymbol{p},\lambda;j})=4i\omega_{p;j}v_{\boldsymbol{p},\lambda;j}^{\dagger}v_{\boldsymbol{p},\lambda;j}=4i\omega_{p;j}\
,$ (38)

where we have made use of the normalization
$v_{\boldsymbol{p},\lambda;j}^{\dagger}v_{\boldsymbol{p},\lambda;j}=1$. Then

$\mathbb{T}^{(N)}_{00}=-2\sum_{\lambda}\sum_{j=1,2}\int d^{3}p\ \omega_{p;j}\
\ \ \ \ \ \ (<0)\ .$ (39)

We can see that $\mathbb{T}^{(N)}_{00}$ corresponds to the diverging negative
energy that is removed by the normal ordering prescription in flat space. In
order that the energy-momentum tensor $T_{\mu\nu}$ approaches the normal-
ordered flat-space expression in the limit, one must subtract
$\mathbb{T}^{(N)}_{\mu\nu}$. Then we define a renormalized energy-momentum
tensor by

$T_{\mu\nu}^{r}=T_{\mu\nu}-\mathbb{T}^{(N)}_{\mu\nu}$ (40)

and whose expectation value is

$\bra{0_{F}(\tau_{0})}T_{\mu\nu}^{r}\ket{0_{F}(\tau_{0})}=\mathbb{T}_{\mu\nu}^{(MIX)}\
.$ (41)

We now derive the general properties of the Bogoliubov coefficients, we show
that $\mathbb{T}^{(MIX)}_{\mu\nu}$ behaves as a perfect fluid and determine
the functional expression of $\mathbb{T}_{\tau\tau}$ and
$\mathbb{T}^{\mu}_{\mu}$.

###  IV.1 General properties of the Bogoliubov Coefficients

Much can still be said without specifying the precise metric, i.e., without
providing an explicit form for the scale factor $C$. Plugging the solutions
(8) in the definition (31) we have

$\displaystyle\Lambda_{p}(\tau)$ $\displaystyle=$ $\displaystyle(2\pi
C)^{3}\left(f_{p;2}^{*}f_{p;1}+g^{*}_{p;2}g_{p;1}\right)\ ,$ (42)

where we have switched to the notation $\Lambda_{p}$ to highlight that it
depends only upon the modulus $p$. In a similar way

$\displaystyle\Xi_{p}(\tau)$ $\displaystyle=$ $\displaystyle(2\pi
C)^{3}\left(f_{p;1}^{*}g^{*}_{p;2}-g^{*}_{p;1}f^{*}_{p;2}\right)\ .$ (43)

From Eqs.(42) and (43) we can verify that

$\displaystyle|\Lambda_{p}(\tau)|^{2}+|\Xi_{p}(\tau)|^{2}$ $\displaystyle=$
$\displaystyle 1\ ,$ (44)

where in the last step we have made use of the normalization condition (72).
We conclude this subsection by expressing the Bogoliubov coefficients in terms
of the reduced functions $\phi_{p},\gamma_{p}$:

$\displaystyle\Lambda_{p}(\tau)$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\left(\phi_{p;2}^{*}\phi_{p;1}+\gamma_{p;2}^{*}\gamma_{p;1}\right)$
$\displaystyle\Xi_{p}(\tau)$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\left(\phi_{p;1}^{*}\gamma^{*}_{p;2}-\gamma_{p;1}^{*}\phi^{*}_{p;2}\right)\
.$ (45)

###  IV.2 Diagonality of the energy-momentum tensor

Using the result (Eq.(IV.1)), we can prove that $\mathbb{T}_{\mu\nu}$ is non-
zero only when $\mu=\nu$. This result has the important consequence that
$\mathbb{T}^{\mu}_{\nu}$ can be interpreted as the energy-momentum tensor of a
perfect fluid.

We start by proving that $\mathbb{T}_{\tau i}=0$ for $i=1,2,3$. Since we can
always write $L_{\tau i}(a,b)=p_{i}h_{a,b}(p)\ ,$ (see Appendix C for
details), each of the terms in the integrand of Eq. (33) is of the form
$p_{i}\mathcal{F}(p)$ , with $\mathcal{F}(p)$ a function of $p$ only. We
therefore have an odd function of $p_{i}$ integrated over an even domain
$p_{i}\in(-\infty,+\infty)$ and the integral vanishes:

$\mathbb{T}_{\tau i}=0\ .$ (46)

The situation is similar for $\mathbb{T}_{ij}$ with $i\neq j$ and $i,j=1,2,3$.
In this case (see Eq.(95) of appendix C), namely,
$L_{ij}(a,b)=p_{i}p_{j}l_{a,b}(p)$, implies that each term under the integral
sign is of the form $p_{i}p_{j}\mathcal{G}(p)$ with $\mathcal{G}(p)$ a
function of $p$ alone. It is clear that also in this case, the integral over
the even domain $(p_{i},p_{j})\in(-\infty,+\infty)\times(-\infty,+\infty)$
vanishes, yielding $\mathbb{T}_{ij}=0$ for $i\neq j$. Notice that the same
conclusion does not apply to $\mathbb{T}_{ii}$, since in this case the
integrand is an even function $p_{i}^{2}\mathcal{A}(p)$ of $p_{i}$. It is easy
to verify that $\mathbb{T}_{ii}$ is the same for each $i=1,2,3$,consistently
with the isotropy of the underlying metric. Due to the manifest diagonality
and isotropy of the energy-momentum tensor, there are really only two
independent components of $\mathbb{T}_{\mu\nu}$, i.e. $\mathbb{T}_{\tau\tau}$
and $\mathbb{T}_{ii}$ for a given $i$. Alternatively one can consider the
$\mathbb{T}_{\tau\tau}$ component and the trace $\mathbb{T}^{\mu}_{\mu}$, as
by definition

$\displaystyle\mathbb{T}_{\mu}^{\mu}$ $\displaystyle=$ $\displaystyle
g^{\mu\nu}\mathbb{T}_{\mu\nu}=C^{-2}\mathbb{T}_{\tau\tau}-3C^{-2}\mathbb{T}_{ii}$

and then

$\displaystyle\mathbb{T}_{ii}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{T}_{\tau\tau}-C^{2}\mathbb{T}^{\mu}_{\mu}}{3}\ ,$
(47)

where no sum is intended over the index $i$ and we have used the isotropy of
$\mathbb{T}_{\mu\nu}$.

It is worth stressing that, as evident from the definition and the expression
of the auxiliary tensor, $\mathbb{T}_{\mu\nu}$ depends only on the time
coordinate $\tau$ and parametrically on the arbitrary fixed time $\tau_{0}$.

We conclude the section by giving the explicit functional form of
$\mathbb{T}_{\tau\tau}$ and $\mathbb{T}^{\mu}_{\mu}$. From Eqs.(33) and (IV.1)
we have

$\displaystyle\mathbb{T}_{\tau\tau}\left[\phi_{p;j},\gamma_{p;j}\right]$
$\displaystyle=$ $\displaystyle 2iC^{-2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\\{|\Xi_{p}(\tau_{0})|^{2}$ (48) $\displaystyle\times$
$\displaystyle\sum_{j=1,2}\left(\phi^{*}_{p;j}\partial_{\tau}\phi_{p;j}+\gamma^{*}_{p;j}\partial_{\tau}\gamma_{p;j}-\partial_{\tau}\phi^{*}_{p;j}\phi_{p;j}-\partial_{\tau}\gamma^{*}_{p;j}\gamma_{p;j}\right)$
$\displaystyle+$ $\displaystyle
2i\imaginary[\Xi_{p}^{*}(\tau_{0})\Lambda_{p}(\tau_{0})\left(\phi^{*}_{p;1}\partial_{\tau}\gamma^{*}_{p;1}-\gamma_{p;1}^{*}\partial_{\tau}\phi^{*}_{p;1}\right)$
$\displaystyle-$
$\displaystyle\Xi^{*}_{p}(\tau_{0})\Lambda^{*}_{p}(\tau_{0})\left(\phi^{*}_{p;2}\partial_{\tau}\gamma^{*}_{p;2}-\gamma_{p;2}^{*}\partial_{\tau}\phi^{*}_{p;2}\right)]\\}$
$\displaystyle-$ $\displaystyle iC^{-2}\sum_{\lambda}\sum_{j=1,2}\int
d^{3}p\left[\phi^{*}_{p;j}\partial_{\tau}\phi_{p;j}+\gamma^{*}_{p;j}\partial_{\tau}\gamma_{p;j}-\partial_{\tau}\phi^{*}_{p;j}\phi_{p;j}-\partial_{\tau}\gamma^{*}_{p;j}\gamma_{p;j}\right]\,,$

and

$\displaystyle\mathbb{T}_{\mu}^{\mu}\left[\phi_{p;j},\gamma_{p;j}\right]$
$\displaystyle=$ $\displaystyle 4iC^{-3}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\\{-i|\Xi_{p}(\tau_{0})|^{2}\sum_{j=1,2}m_{j}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)$
(49) $\displaystyle+$ $\displaystyle
2i\imaginary\left[im_{2}\Xi_{p}^{*}(\tau_{0})\Lambda^{*}_{p}(\tau_{0})\phi^{*}_{p;2}\gamma^{*}_{p;2}-im_{1}\Xi_{p}^{*}(\tau_{0})\Lambda_{p}(\tau_{0})\phi^{*}_{p;1}\gamma^{*}_{p;1}\right]\\}$
$\displaystyle-$ $\displaystyle 2C^{-3}\sum_{\lambda}\sum_{j=1,2}m_{j}\int
d^{3}p\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)\ .$

Note that both the $\tau\tau$ component (48) and the trace (49) can also be
seen as functionals of $\\{\phi_{p},\partial_{\tau}\phi_{p}\\}$. In fact from
Eqs.(15) we have

$\displaystyle\gamma_{p}(\eta)$ $\displaystyle=$
$\displaystyle\frac{i}{p}\left(\partial_{\eta}\phi_{p}(\eta)+imC(\eta)\phi_{p}(\eta)\right)$
$\displaystyle\partial_{\eta}\gamma_{p}(\eta)$ $\displaystyle=$
$\displaystyle\frac{i}{p}\left(\partial_{\eta}^{2}\phi_{p}(\eta)+im\partial_{\eta}C(\eta)\phi_{p}(\eta)+imC(\eta)\partial_{\eta}\phi_{p}(\eta)\right)\
,$ (50)

which can be substituted in the expressions above for $\mathbb{T}_{\mu\nu}$.
From the diagonality of the expectation value it is also straightforward to
prove its covariant conservation $\nabla_{\mu}\mathbb{T}^{\mu\nu}=0$ (see the
appendix E for the details).

##  V Applications: exponential expansion

We have shown that independently of the specific scale factor $C$, the energy-
momentum tensor associated to the flavor vacuum $\mathbb{T}_{\mu\nu}^{(MIX)}$
satisfies a number of important properties: (i) it is diagonal, (ii) it is
covariantly conserved and (iii) depends only on time $\tau$. Then
$\mathbb{T}_{\mu\nu}^{(MIX)}$ for the metric (1) corresponds to the energy-
momentum tensor of a perfect fluid with time-dependent energy density and
pressure. In this section, in order to better understand the properties of
$\mathbb{T}_{\mu\nu}^{(MIX)}$, we assume a specific evolution of the scale
factor $C$ and compute the corresponding expectation value of the energy-
momentum tensor on the flavor vacuum. In doing so, we are neglecting the back-
reaction due to the flavor fields, i.e., the modifications induced on the
metric by the energy-momentum tensor of Eq. (27). The computation based on a
fixed background metric is undoubtedly an approximation, but it is useful to
get an insight into the kind of contribution that emerges from the flavor
vacuum. The self-consistent way to deal with $\mathbb{T}_{\mu\nu}$ is to
insert Eq. (27), together with all the relevant matter terms, on the right-
hand side of the Einstein equations and then solve simultaneously for the
scale factor $C$ and the Dirac modes. This kind of calculation will be
performed elsewhere.

###  V.1 Positive energy solutions

Here, in particular, we study the energy momentum-tensor corresponding to an
exponential evolution of the scale factor $C(t)=e^{H_{0}t}$, with $H_{0}$ a
constant with dimensions of mass. The great advantage of so-picked scale
factor is that the corresponding Dirac equation can be solved analytically
without resorting to any approximation, allowing for an in-depth analysis of
$\mathbb{T}^{(MIX)}_{\mu\nu}$. Transforming to conformal time we have

$\tau=-\frac{1}{H_{0}}e^{-H_{0}t}\,,\ \ \ \ \ \ \ \ \ C=-\frac{1}{H_{0}\tau}\
.$ (51)

Notice that the conformal time $\tau$ is always negative. Inserting Eq. (51)
in Eq. (17) we obtain

$\tau^{2}\partial_{\tau}^{2}\phi_{p}+\left(p^{2}\tau^{2}+\frac{im}{H_{0}}+\frac{m^{2}}{H_{0}^{2}}\right)\phi_{p}=0\
,$ (52)

or, introducing the positive variable $s=-p\tau$

$s^{2}\partial_{s}^{2}\phi_{p}+\left(s^{2}+\frac{im}{H_{0}}+\frac{m^{2}}{H_{0}^{2}}\right)\phi_{p}=0\
.$ (53)

This is a Bessel-like equation, whose general solution can be written as

$\phi_{p}(s)=s^{\frac{1}{2}}\left(C_{1}J_{\nu}(s)+C_{2}J_{-\nu}(s)\right)\,,\
\ \ \ \ \ \ \ \ \ \nu=\frac{1}{2}\left(1-\frac{2im}{H_{0}}\right)$ (54)

with $J_{\nu}(s)$ denoting the Bessel function of order $\nu$ and
$C_{1},C_{2}$ arbitrary complex constants. In order to specify the solution we
need to impose some kind of boundary conditions, which in turn determine the
positive energy solutions. We require that the modes of Eq. (54) be positive
energy with respect to $\partial_{\tau}$ at early times, i.e. for
$\tau\rightarrow-\infty$ (where $C\rightarrow 0$). With this choice the mass
vacuum corresponds to the absence of massive neutrinos at early times. As
$\tau\rightarrow-\infty$ the mass terms can be neglected, and equation (52)
becomes

$\partial^{2}_{\tau}\phi_{p}+p^{2}\phi_{p}=0\ .$ (55)

The positive energy solution with respect to $\partial_{\tau}$ is evidently
$\phi^{+}_{p}\propto e^{-ip\tau}$, or, in terms of the $s$ variable,
$\phi^{+}_{p}(s)\propto e^{is}$. We then impose the requirement
$\lim_{s\rightarrow+\infty}\phi_{p}(s)\propto e^{is}\ .$ Recalling that $s$ is
a positive real variable, we can employ the large argument expansion of the
Bessel functions [57]

$J_{\nu}(s)\simeq\sqrt{\frac{2}{\pi
s}}\cos\left(s-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)\ \ \ \ for\ \ \ \ \
s\rightarrow+\infty$ (56)

In this way, we show that the combination satisfying the requirement is given
by

$\phi_{p}(s)=N_{p}s^{\frac{1}{2}}\left(J_{\nu}(s)-ie^{\frac{\pi
m}{H_{0}}}J_{-\nu}(s)\right)\simeq N_{p}\sqrt{\frac{2}{\pi}}(-ie^{\frac{\pi
m}{2H_{0}}})\cosh\left(\frac{\pi m}{H_{0}}\right)e^{is}\ ,$ (57)

where $N_{p}$ is a normalization constant and the last equivalence holds in
the limit $s\rightarrow+\infty$. Inserting Eq.(57) in the first of the
Eqs.(IV.2), we deduce

$\gamma_{p}(s)=N_{p}s^{\frac{1}{2}}\left(-iJ_{\nu-1}(s)+e^{\frac{\pi
m}{H_{0}}}J_{1-\nu}(s)\right)$ (58)

where we have made use of the differential relations satisfied by the Bessel
functions [57]. In order to fix $N_{p}$ we impose the normalization condition
(16) in the $s\rightarrow+\infty$ limit. Using the large argument expansion
once again, we obtain $|N_{p}|^{2}=\frac{1}{32\pi^{2}\cosh^{2}\left(\frac{\pi
m}{H_{0}}\right)}e^{-\frac{\pi m}{H_{0}}}$. Finally, choosing $N_{p}$ real, we
can write the positive energy solutions as

$\displaystyle\phi_{p}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\frac{e^{-\frac{\pi
m}{2H_{0}}}}{\cosh\left(\frac{\pi
m}{H_{0}}\right)}\sqrt{\frac{s}{2}}\left(J_{\nu}(s)-ie^{\frac{\pi
m}{H_{0}}}J_{-\nu}(s)\right)$ $\displaystyle\gamma_{p}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\frac{e^{-\frac{\pi
m}{2H_{0}}}}{\cosh\left(\frac{\pi
m}{H_{0}}\right)}\sqrt{\frac{s}{2}}\left(-iJ_{\nu-1}(s)+e^{\frac{\pi
m}{H_{0}}}J_{1-\nu}(s)\right)$ (59)

with $\nu=\frac{1}{2}\left(1-\frac{2im}{H_{0}}\right)$ and $s=-p\tau$.

###  V.2 Bogoliubov Coefficients

We have two sets of solutions $\phi_{p,j},\gamma_{p,j}$, one for each value of
the mass $m_{j}$, with $j=1,2$ from the equations (V.1). The compatibility
requirement then implies that each of the $\phi_{p,j},\gamma_{p,j}$ has the
same form for $j=1,2$, except that one has $m_{j}$ wherever the mass appears,
including the function index
$\nu_{j}=\frac{1}{2}\left(1-\frac{2im_{j}}{H_{0}}\right)$. We can compute the
Bogoliubov coefficients ($\Lambda_{p}(\tau),\Xi_{p}(\tau)$) straight away from
Eqs.(IV.1):

$\displaystyle\Lambda_{p}(\tau)$ $\displaystyle=$ $\displaystyle\frac{\pi
s}{4}\frac{e^{-\frac{\pi}{2H_{0}}(m_{1}+m_{2})}}{\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\Bigg{\\{}J_{\nu_{2}}^{*}(s)J_{\nu_{1}}(s)+J^{*}_{\nu_{2}-1}(s)J_{\nu_{1}-1}(s)+ie^{\frac{\pi
m_{1}}{H_{0}}}\left[J^{*}_{\nu_{2}-1}(s)J_{1-\nu_{1}}(s)-J^{*}_{\nu_{2}}(s)J_{-\nu_{1}}(s)\right]$
$\displaystyle+$ $\displaystyle ie^{\frac{\pi
m_{2}}{H_{0}}}\left[J^{*}_{-\nu_{2}}(s)J_{\nu_{1}}(s)-J^{*}_{1-\nu_{2}}(s)J_{\nu_{1}-1}(s)\right]+e^{\frac{\pi}{H_{0}}(m_{1}+m_{2})}\left[J^{*}_{-\nu_{2}}(s)J_{-\nu_{1}}(s)+J^{*}_{1-\nu_{2}}(s)J_{1-\nu_{1}}(s)\right]\Bigg{\\}}$
$\displaystyle\Xi_{p}(\tau)$ $\displaystyle=$ $\displaystyle\frac{\pi
s}{4}\frac{e^{-\frac{\pi}{2H_{0}}(m_{1}+m_{2})}}{\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\Bigg{\\{}i\left[J^{*}_{\nu_{1}}(s)J^{*}_{\nu_{2}-1}(s)-J^{*}_{\nu_{1}-1}(s)J^{*}_{\nu_{2}}(s)\right]+e^{\frac{\pi
m_{2}}{H_{0}}}\left[J^{*}_{\nu_{1}}(s)J^{*}_{1-\nu_{2}}(s)+J^{*}_{\nu_{1}-1}(s)J^{*}_{-\nu_{2}}(s)\right]$
(60) $\displaystyle-$ $\displaystyle e^{\frac{\pi
m_{1}}{H_{0}}}\left[J^{*}_{-\nu_{1}}(s)J^{*}_{\nu_{2}-1}(s)+J^{*}_{1-\nu_{1}}(s)J^{*}_{\nu_{2}}(s)\right]+ie^{\frac{\pi}{H_{0}}(m_{1}+m_{2})}\left[J^{*}_{-\nu_{1}}(s)J^{*}_{1-\nu_{2}}(s)-J^{*}_{1-\nu_{1}}(s)J^{*}_{-\nu_{2}}(s)\right]\Bigg{\\}}$

To give a flavor of the behaviour of the Bogoliubov coefficients we have
plotted $|\Xi_{p}(s)|^{2}$ as a function of $s$ for sample values of
$m_{1},m_{2}$ in figure (1).

Figure 1: Squared modulus of the Bogoliubov coefficient $|\Xi_{p}|^{2}$ as a
function of $s$ for sample values of the masses (all in units of $H_{0}$):
(black dotted line) $m_{1}=0.7$, $m_{2}=1.4$ , (red dashed line)
$m_{1}=1,m_{2}=2$, (blue dotdashed line) $m_{1}=10,m_{2}=20$, (dark orange
solid line), $m_{1}=100$, $m_{2}=300$. The momentum dependence is implicit in
$s=-p\tau$.

We shall be particularly interested in the late time expression of the
Bogoliubov coefficients, namely for $s\rightarrow 0^{+}$ (the corresponding
limit is $\tau\rightarrow 0^{-}$, or $t\rightarrow+\infty$). For its
determination we make use of the small argument expansion of the Bessel
functions
$J_{\nu}(s)\simeq\left(\frac{s}{2}\right)^{\nu}\frac{1}{\Gamma(1+\nu)}$ with
$\Gamma$ denoting the Euler gamma function. From the Eqs.(V.2), it is easy to
find that at the leading order for $s\rightarrow 0^{+}$ one has

$\displaystyle\Lambda_{p}(s)$ $\displaystyle\simeq$
$\displaystyle\frac{\pi}{2}\frac{e^{-\frac{\pi}{2H_{0}}(m_{1}+m_{2})}}{\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\left[\frac{e^{i\left(\frac{m_{2}-m_{1}}{H_{0}}\right)\log\left(\frac{s}{2}\right)}}{\Gamma^{*}(\nu_{2})\Gamma(\nu_{1})}+e^{\frac{\pi}{H_{0}}(m_{1}+m_{2})}\frac{e^{-i\left(\frac{m_{2}-m_{1}}{H_{0}}\right)\log\left(\frac{s}{2}\right)}}{\Gamma^{*}(1-\nu_{2})\Gamma(1-\nu_{1})}\right]$
$\displaystyle\Xi_{p}(s)$ $\displaystyle\simeq$
$\displaystyle\frac{\pi}{2}\frac{e^{-\frac{\pi}{2H_{0}}(m_{1}+m_{2})}}{\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\left[\frac{e^{\frac{\pi
m_{2}}{H_{0}}}e^{-i\left(\frac{m_{2}-m_{1}}{H_{0}}\right)\log\left(\frac{s}{2}\right)}}{\Gamma^{*}(\nu_{1})\Gamma^{*}(1-\nu_{2})}-\frac{e^{\frac{\pi
m_{1}}{H_{0}}}e^{i\left(\frac{m_{2}-m_{1}}{H_{0}}\right)\log\left(\frac{s}{2}\right)}}{\Gamma^{*}(1-\nu_{1})\Gamma^{*}(\nu_{2})}\right]$
(61)

###  V.3 Explicit form of the energy-momentum tensor

For the explicit calculation it is convenient to refer to the splitting of
Eq.(33) and compute $\mathbb{T}_{\mu\nu}^{(MIX)}$ and
$\mathbb{T}_{\mu\nu}^{(N)}$ separately. Moreover, keeping in mind the results
of the previous sections, it is sufficient to compute the $\tau\tau$ component
and the trace in order to fully determine $\mathbb{T}_{\mu\nu}$. Inserting the
solutions (V.1) in the Eq.(48) we find, after a lengthy but straightforward
calculation:

$\displaystyle\mathbb{T}^{(MIX)}_{\tau\tau}(\tau)$ $\displaystyle=$
$\displaystyle i\sin^{2}\theta\sum_{\lambda}\int
d^{3}p|\Xi_{p}(\tau_{0})|^{2}\left(\frac{H^{2}_{0}p^{2}\tau^{3}}{16\pi^{2}}\right)\sum_{j=1,2}\frac{e^{-\frac{\pi
m_{j}}{H_{0}}}}{\cosh^{2}\left(\frac{\pi
m_{j}}{H_{0}}\right)}\Bigg{\\{}\bigg{[}2(J^{*}_{\nu_{j}}J_{\nu_{j}-1}-J^{*}_{\nu_{j}-1}J_{\nu})+\frac{\nu_{j}^{*}-\nu_{j}}{s}|J_{\nu_{j}}|^{2}$
$\displaystyle+$
$\displaystyle\frac{\nu_{j}-\nu_{j}^{*}}{s}|J_{\nu_{j}-1}|^{2}\bigg{]}+ie^{\frac{\pi
m_{j}}{H_{0}}}\bigg{[}2\left(J^{*}_{\nu_{j}}J_{1-\nu_{j}}+J^{*}_{-\nu_{j}}J_{\nu_{j}-1}+J^{*}_{\nu_{j}-1}J_{-\nu_{j}}+J^{*}_{1-\nu_{j}}J_{\nu_{j}}\right)+\frac{\nu_{j}-\nu_{j}^{*}}{s}J^{*}_{\nu_{j}}J_{-\nu_{j}}$
$\displaystyle+$
$\displaystyle\frac{\nu_{j}^{*}-\nu_{j}}{s}J^{*}_{-\nu_{j}}J_{\nu_{j}}+\frac{\nu_{j}-\nu_{j}^{*}}{s}J^{*}_{\nu_{j}-1}J_{1-\nu_{j}}+\frac{\nu_{j}^{*}-\nu_{j}}{s}J^{*}_{1-\nu_{j}}J_{\nu_{j}-1}\bigg{]}+e^{\frac{2\pi
m_{j}}{H_{0}}}\bigg{[}2(J^{*}_{1-\nu_{j}}J_{-\nu_{j}}-J^{*}_{-\nu_{j}}J_{1-\nu_{j}})$
$\displaystyle+$
$\displaystyle\frac{\nu_{j}^{*}-\nu_{j}}{s}|J_{-\nu_{j}}|^{2}+\frac{\nu_{j}-\nu_{j}^{*}}{s}|J_{1-\nu_{j}}|^{2}\bigg{]}\Bigg{\\}}$
$\displaystyle+$ $\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\Bigg{\\{}\Xi^{*}_{p}(\tau_{0})\Lambda_{p}(\tau_{0})\left(\frac{H_{0}^{2}p^{2}\tau^{3}e^{-\frac{\pi
m_{1}}{H_{0}}}}{8\pi^{2}\cosh^{2}\left(\frac{\pi
m_{1}}{H_{0}}\right)}\right)\bigg{[}\left(-i(J^{*}_{\nu_{1}})^{2}-i(J^{*}_{\nu_{1}-1})^{2}+i\frac{2\nu^{*}_{1}-1}{s}J^{*}_{\nu_{1}}J^{*}_{\nu_{1}-1}\right)$
$\displaystyle+$ $\displaystyle e^{\frac{\pi
m_{1}}{H_{0}}}\left(2(J^{*}_{\nu_{1}}J^{*}_{-\nu_{1}}-J^{*}_{\nu_{1}-1}J^{*}_{1-\nu_{1}})+\frac{2\nu_{1}^{*}-1}{s}J^{*}_{\nu_{1}}J^{*}_{1-\nu_{1}}+\frac{1-2\nu_{1}^{*}}{s}J^{*}_{-\nu_{1}}J^{*}_{\nu_{1}-1}\right)$
$\displaystyle+$ $\displaystyle ie^{\frac{2\pi
m_{1}}{H_{0}}}\left((J^{*}_{-\nu_{1}})^{2}+(J^{*}_{1-\nu_{1}})^{2}+\frac{2\nu_{1}^{*}-1}{s}J^{*}_{-\nu_{1}}J^{*}_{1-\nu_{1}}\right)\bigg{]}-c.c.\Bigg{\\}}$
$\displaystyle-$ $\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\Bigg{\\{}\Xi^{*}_{p}(\tau_{0})\Lambda^{*}_{p}(\tau_{0})\left(\frac{H_{0}^{2}p^{2}\tau^{3}e^{-\frac{\pi
m_{2}}{H_{0}}}}{8\pi^{2}\cosh^{2}\left(\frac{\pi
m_{2}}{H_{0}}\right)}\right)\bigg{[}\left(-i(J^{*}_{\nu_{2}})^{2}-i(J^{*}_{\nu_{2}-1})^{2}+i\frac{2\nu^{*}_{2}-1}{s}J^{*}_{\nu_{2}}J^{*}_{\nu_{2}-1}\right)$
$\displaystyle+$ $\displaystyle e^{\frac{\pi
m_{2}}{H_{0}}}\left(2(J^{*}_{\nu_{2}}J^{*}_{-\nu_{2}}-J^{*}_{\nu_{2}-1}J^{*}_{1-\nu_{2}})+\frac{2\nu_{2}^{*}-1}{s}J^{*}_{\nu_{2}}J^{*}_{1-\nu_{2}}+\frac{1-2\nu_{2}^{*}}{s}J^{*}_{-\nu_{2}}J^{*}_{\nu_{2}-1}\right)$
$\displaystyle+$ $\displaystyle ie^{\frac{2\pi
m_{2}}{H_{0}}}\left((J^{*}_{-\nu_{2}})^{2}+(J^{*}_{1-\nu_{2}})^{2}+\frac{2\nu_{2}^{*}-1}{s}J^{*}_{-\nu_{2}}J^{*}_{1-\nu_{2}}\right)\bigg{]}-c.c.\Bigg{\\}}$
(62)

In the above equation we have left the Bogoliubov coefficients implicit and
the argument of the Bessel functions $s=-p\tau$ has been omitted for a better
visualization. We recall that the helicity sum is over $\lambda=\pm 1$ and
that in general $\tau\neq\tau_{0}$. An analogous calculation can be performed
for the trace:

$\displaystyle\mathbb{T}^{\mu(MIX)}_{\mu}$ $\displaystyle=$ $\displaystyle
i\sin^{2}\theta\sum_{\lambda}\int
d^{3}p|\Xi_{p}(\tau_{0})|^{2}\left(\frac{iH_{0}^{3}\tau^{3}s}{8\pi^{2}}\right)\sum_{j=1,2}\left(\frac{m_{j}e^{\frac{-\pi
m_{j}}{H_{0}}}}{\cosh^{2}\left(\frac{\pi
m_{j}}{H_{0}}\right)}\right)\Bigg{\\{}|J_{\nu_{j}}|^{2}-|J_{\nu_{j}-1}|^{2}$
(63) $\displaystyle+$ $\displaystyle ie^{\frac{\pi
m_{j}}{H_{0}}}\left(J^{*}_{-\nu_{j}}J_{\nu_{j}}-J^{*}_{\nu_{j}}J_{-\nu_{j}}+J^{*}_{1-\nu_{j}}J_{\nu_{j}-1}-J^{*}_{\nu_{j}-1}J_{1-\nu_{j}}\right)+e^{\frac{2\pi
m_{j}}{H_{0}}}\left(|J_{-\nu_{j}}|^{2}-|J_{1-\nu_{j}}|^{2}\right)\Bigg{\\}}$
$\displaystyle+$ $\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\Bigg{\\{}\Xi_{p}^{*}(\tau_{0})\Lambda_{p}(\tau_{0})\left(\frac{im_{1}sH_{0}^{3}\tau^{3}e^{-\frac{\pi
m_{1}}{H_{0}}}}{4\pi^{2}\cosh^{2}\left(\frac{\pi
m_{1}}{H_{0}}\right)}\right)\bigg{[}iJ^{*}_{\nu_{1}}J^{*}_{\nu_{1}-1}+e^{\frac{\pi
m_{1}}{H_{0}}}J^{*}_{\nu_{1}}J^{*}_{1-\nu_{1}}$ $\displaystyle-$
$\displaystyle e^{\frac{\pi
m_{1}}{H_{0}}}J^{*}_{-\nu_{1}}J^{*}_{\nu_{1}-1}+ie^{\frac{2\pi
m_{1}}{H_{0}}}J^{*}_{-\nu_{1}}J^{*}_{1-\nu_{1}}\bigg{]}-c.c\Bigg{\\}}$
$\displaystyle-$ $\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\Bigg{\\{}\Xi_{p}^{*}(\tau_{0})\Lambda^{*}_{p}(\tau_{0})\left(\frac{im_{2}sH_{0}^{3}\tau^{3}e^{-\frac{\pi
m_{2}}{H_{0}}}}{4\pi^{2}\cosh^{2}\left(\frac{\pi
m_{2}}{H_{0}}\right)}\right)\bigg{[}iJ^{*}_{\nu_{2}}J^{*}_{\nu_{2}-1}+e^{\frac{\pi
m_{2}}{H_{0}}}J^{*}_{\nu_{2}}J^{*}_{1-\nu_{2}}$ $\displaystyle-$
$\displaystyle e^{\frac{\pi
m_{2}}{H_{0}}}J^{*}_{-\nu_{2}}J^{*}_{\nu_{2}-1}+ie^{\frac{2\pi
m_{2}}{H_{0}}}J^{*}_{-\nu_{2}}J^{*}_{1-\nu_{2}}\bigg{]}-c.c\Bigg{\\}}\ \ \ \ \
\ \ .$

We are particularly interested in the late time ($\tau\rightarrow 0^{-}$)
expression of the above equations. According to the definition (27), these
represent the contribution of the flavor vacuum state, defined at an earlier
time $\tau_{0}<\tau$, to the energy-momentum tensor at late times. We then
perform the small argument expansion
$J_{\nu}(-p\tau)\simeq\left(\frac{-p\tau}{2}\right)^{\nu}\frac{1}{\Gamma(1+\nu)}$
for all the Bessel functions appearing in Eqs. (V.3) and (63), and keep only
the terms of lowest order in the variable $\tau$. At order $\tau$, the
$\tau\tau$ component is found to be

$\displaystyle\mathbb{T}_{\tau\tau}^{(MIX)(1)}(\tau)$ $\displaystyle\simeq$
$\displaystyle i\sin^{2}\theta\sum_{\lambda}\int
d^{3}p|\Xi_{p}(\tau_{0})|^{2}\left(i\frac{H_{0}\tau}{2\pi^{3}}\right)\sum_{j=1,2}m_{j}\tanh\left(\frac{\pi
m_{j}}{H_{0}}\right)$ (64) $\displaystyle+$
$\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\left[\Xi^{*}_{p}(\tau_{0})\Lambda_{p}(\tau_{0})\left(\frac{-im_{1}H_{0}\tau}{2\pi^{3}\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)}\right)-c.c.\right]$ $\displaystyle-$
$\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\left[\Xi^{*}_{p}(\tau_{0})\Lambda^{*}_{p}(\tau_{0})\left(\frac{-im_{2}H_{0}\tau}{2\pi^{3}\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\right)-c.c.\right]\ .$

The corresponding lowest order in the trace is $\propto\tau^{3}$. To see why
this is the case, recall that by definition

$\mathbb{T}_{\mu}^{\mu}=C^{-2}\mathbb{T}_{\tau\tau}-C^{-2}\sum_{l=1}^{3}\mathbb{T}_{ll}=H_{0}^{2}\tau^{2}\mathbb{T}_{\tau\tau}-H_{0}^{2}\tau^{2}\sum_{l=1}^{3}\mathbb{T}_{ll}\
,$ (65)

so that in correspondence with $\mathbb{T}_{\tau\tau}\propto\tau$ one has
$\mathbb{T}^{\mu}_{\mu}\propto\tau^{3}$. To this order the trace is

$\displaystyle\mathbb{T}^{\mu\ (MIX)\ (1)}_{\mu}$ $\displaystyle\simeq$
$\displaystyle i\sin^{2}\theta\sum_{\lambda}\int
d^{3}p|\Xi_{p}(\tau_{0})|^{2}\left(\frac{iH_{0}^{3}\tau^{3}}{2\pi^{3}}\right)\sum_{j=1,2}m_{j}\tanh\left(\frac{\pi
m_{j}}{H_{0}}\right)$ (66) $\displaystyle+$
$\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\left[\Xi^{*}_{p}(\tau_{0})\Lambda_{p}(\tau_{0})\left(\frac{-im_{1}H_{0}^{3}\tau^{3}}{2\pi^{3}\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)}\right)-c.c.\right]$ $\displaystyle-$
$\displaystyle\frac{i}{2}\sin^{2}\theta\sum_{\lambda}\int
d^{3}p\left[\Xi^{*}_{p}(\tau_{0})\Lambda^{*}_{p}(\tau_{0})\left(\frac{-im_{2}H_{0}^{3}\tau^{3}}{2\pi^{3}\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\right)-c.c.\right]\ .$

Inserting Eqs. (64) and (66) in Eq. (47), we can deduce the important result
$\mathbb{T}_{ii}^{(MIX)(1)}=0\ \,,\forall i\ .$ In other words, at first order
in $\tau$, the spatial components of the energy-momentum tensor vanish. Then
the equation of state reads, at lowest order in $\tau$:

$w^{(1)}=\frac{\mathbb{T}^{i\ (MIX)(1)}_{i}}{\mathbb{T}^{\tau\
(MIX)(1)}_{\tau}}=0\ ,$ (67)

i. e., _the energy-momentum tensor associated to the flavor vacuum satisfies,
at late times ( $\tau\rightarrow 0^{-}$), the equation of state of a pressure-
less perfect fluid_. It is important to stress that this result does not
depend on the value of the momentum integrals, and therefore is independent of
any regularization. Moreover, as it is evident from (67) it holds for any
choice of the reference time $\tau_{0}$.

###  V.4 Regularization

All the momentum integrals in the expression for the energy-momentum tensor,
both for the general case and the for the late time approximation, are to be
performed over the whole of $\mathbb{R}^{3}$ and are formally divergent. To
extract a finite result, we need some form of regularization. The most
immediate way to regularize the momentum integrations is to introduce an
ultraviolet momentum cut-off $\mathcal{K}$. Usually, the cut-off is chosen in
correspondence with a ’new physics’ energy scale, beyond which the ’low
energy’ quantum field theory description breaks down. In our case,
$\mathcal{K}$ ought to be related to the scale at which the semiclassical
approximation breaks down, i. e., at the quantum gravity scale
$\mathcal{K}\simeq M_{pl}$ which is of the order of the Planck mass $M_{pl}$.
However, we can expect the cut-off to be at a much lower scale for what
concerns the proper mixing term. Indeed, at very high energies, the mass
difference between the neutrino states becomes negligible, and the oscillation
frequency $\propto\frac{1}{p}$ approaches zero, implying that there is no
oscillation. The same result also holds in quantum field theory, where the
mixing Bogoliubov coefficient $\Xi_{p}$ generally approaches zero at high
energies (therefore yielding no contribution to the energy-momentum tensor.).
The heuristic argument above provides a physical reason to adopt the cut-off
regularization and also justifies the adoption of a cut-off scale
$\mathcal{K}\ll M_{pl}$, at least for what concerns the mixing.

Before proceeding we need to clarify that the cut-off must be imposed upon the
comoving momentum $\boldsymbol{p}_{PHYS}=\frac{\boldsymbol{p}}{C}\ .$ rather
than the mode label $p$.

The imposition of a cut-off $\mathcal{K}_{0}$ on $p_{PHYS}$ then translates
into a sort of ’comoving’ cutoff for the mode label $p$

$p_{CUTOFF}=\mathcal{K}_{0}C=\frac{-\mathcal{K}_{0}}{H_{0}\tau}=\mathcal{K}(\tau)\
,$ (68)

which is strictly positive (recall that $\tau<0$). For the late time energy-
momentum tensor, in the approximation in which $\tau_{0}<\tau$ is also at late
times $\tau_{0}\rightarrow 0^{-}$, we can give a simple analytical form of the
regularized integrals. Performing the integrals in Eqs. (64) and (66), with
the comoving cut-off of Eq. (68), we obtain

$\displaystyle\mathbb{T}_{\tau\tau}^{(MIX)\ (1)}$ $\displaystyle=$
$\displaystyle\frac{-\sin^{2}\theta
H_{0}\tau\mathcal{K}^{3}(\tau)}{3\pi^{2}}\Bigg{[}\frac{e^{\frac{\pi(m_{2}-m_{1})}{H_{0}}}+e^{-\frac{\pi(m_{2}-m_{1})}{H_{0}}}}{\cosh\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh\left(\frac{\pi
m_{2}}{H_{0}}\right)}\left(m_{1}\tanh\left(\frac{\pi
m_{1}}{H_{0}}\right)+m_{2}\tanh\left(\frac{\pi m_{2}}{H_{0}}\right)\right)$
$\displaystyle-$ $\displaystyle 2\frac{m_{1}\tanh\left(\frac{\pi
m_{2}}{H_{0}}\right)}{\cosh^{2}\left(\frac{\pi
m_{1}}{H_{0}}\right)}-2\frac{m_{2}\tanh\left(\frac{\pi
m_{1}}{H_{0}}\right)}{\cosh^{2}\left(\frac{\pi
m_{2}}{H_{0}}\right)}\Bigg{]}+\sin^{2}\theta
H_{0}\tau\mathcal{K}^{3}(\tau)\Bigg{[}\left(\frac{m_{1}\tanh\left(\frac{\pi
m_{1}}{H_{0}}\right)+m_{2}\tanh\left(\frac{\pi
m_{2}}{H_{0}}\right)}{\cosh^{2}\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh^{2}\left(\frac{\pi m_{2}}{H_{0}}\right)}\right)$
$\displaystyle\times$
$\displaystyle\left(\frac{1}{\Gamma(\nu_{1})\Gamma^{*}(\nu_{2})}\right)^{2}\left(\frac{1}{3+2i\frac{m_{2}-m_{1}}{H_{0}}}\right)\left(\frac{-\mathcal{K}(\tau)\tau_{0}}{2}\right)^{\frac{2i(m_{2}-m_{1})}{H_{0}}}+c.c.\Bigg{]}$
$\displaystyle+$ $\displaystyle
i\sin^{2}\theta\mathcal{K}^{3}(\tau)\Bigg{\\{}\Bigg{[}\left(\frac{-im_{1}H_{0}\tau
e^{\frac{-\pi m_{1}}{H_{0}}}}{2\cosh^{3}\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh^{2}\left(\frac{\pi
m_{2}}{H_{0}}\right)}\right)\left(\frac{1}{\Gamma(\nu_{1})\Gamma^{*}(\nu_{2})}\right)^{2}\left(\frac{1}{3+2i\frac{m_{2}-m_{1}}{H_{0}}}\right)\left(\frac{-\mathcal{K}(\tau)\tau_{0}}{2}\right)^{\frac{2i(m_{2}-m_{1})}{H_{0}}}$
$\displaystyle-$ $\displaystyle\left(\frac{-im_{1}H_{0}\tau e^{\frac{-\pi
m_{1}}{H_{0}}}}{2\cosh^{3}\left(\frac{\pi
m_{1}}{H_{0}}\right)\cosh^{2}\left(\frac{\pi
m_{2}}{H_{0}}\right)}\right)\left(\frac{1}{\Gamma^{*}(\nu_{1})\Gamma(\nu_{2})}\right)^{2}\left(\frac{1}{3-2i\frac{m_{2}-m_{1}}{H_{0}}}\right)\left(\frac{-\mathcal{K}(\tau)\tau_{0}}{2}\right)^{\frac{-2i(m_{2}-m_{1})}{H_{0}}}\Bigg{]}-c.c.\Bigg{\\}}$
$\displaystyle-$ $\displaystyle
i\sin^{2}\theta\mathcal{K}^{3}(\tau)\Bigg{\\{}\Bigg{[}\left(\frac{-im_{2}H_{0}\tau
e^{\frac{-\pi m_{2}}{H_{0}}}}{2\cosh^{3}\left(\frac{\pi
m_{2}}{H_{0}}\right)\cosh^{2}\left(\frac{\pi
m_{1}}{H_{0}}\right)}\right)\left(\frac{1}{\Gamma(\nu_{1})\Gamma^{*}(\nu_{2})}\right)^{2}\left(\frac{1}{3+2i\frac{m_{2}-m_{1}}{H_{0}}}\right)\left(\frac{-\mathcal{K}(\tau)\tau_{0}}{2}\right)^{\frac{2i(m_{2}-m_{1})}{H_{0}}}$
$\displaystyle-$ $\displaystyle\left(\frac{-im_{2}H_{0}\tau e^{\frac{-\pi
m_{2}}{H_{0}}}}{2\cosh^{3}\left(\frac{\pi
m_{2}}{H_{0}}\right)\cosh^{2}\left(\frac{\pi
m_{1}}{H_{0}}\right)}\right)\left(\frac{1}{\Gamma^{*}(\nu_{1})\Gamma(\nu_{2})}\right)^{2}\left(\frac{1}{3-2i\frac{m_{2}-m_{1}}{H_{0}}}\right)\left(\frac{-\mathcal{K}(\tau)\tau_{0}}{2}\right)^{\frac{-2i(m_{2}-m_{1})}{H_{0}}}\Bigg{]}-c.c.\Bigg{\\}}$
$\displaystyle\ \ \ \ \ \ .$ (69)

The trace can be deduced immediately from Eq. (V.4) by simply multiplying by
the factor $C^{-2}$, since the pressure is zero for the late time energy
momentum tensor (see Eq. (67)).

Figure 2:  Logarithmic scale plot of the energy density
$\rho_{MIX}=\mathbb{T}^{\tau(MIX)}_{\tau}$ from Eq. (V.4) as a function of
conformal time $\tau$ for sample values of the parameters. The corresponding
coordinate time $t$ is reported above. We have used a cut-off
$\mathcal{K}_{0}=246\ \mathrm{GeV}$ of the order of the electroweak scale,
neutrino masses $m_{1}=15.25H_{0}\ ,m_{2}=22.25H_{0}$ and the expansion rate
$H_{0}=10^{-3}\mathrm{eV}$.

A visual indication of the physical content of Eq. (V.4) is provided in Figure
(2), where the mixing energy density $\mathbb{T}^{\tau\ (MIX)}_{\tau}$ from
equation (V.4) is plotted against conformal time $\tau$ for sample values of
the parameters. We notice that $\rho_{MIX}$, for the parameters and the range
considered in Figure (2), is nearly constant, except for tiny oscillations of
relative magnitude $\delta\rho_{MIX}/\rho_{MIX}\simeq 10^{-9}$. The
oscillations come from the imaginary exponentials in Eq. (V.4). Such a simple
evolution pattern can be expected to disappear when other time ranges are
considered, and the full expression from Eq. (V.3) is employed, entailing a
much more intricate time evolution. The results obtained in this work indicate
that the vacuum energy associate to neutrino mixing in curved space might
represent a dark matter component. Notice that for a momentum cutoff of the
order of the electroweak scale $\mathcal{K}_{0}=246\ \mathrm{GeV}$, for
neutrino masses $m_{1},m_{2}$ such that $\Delta m^{2}_{12}\simeq
10^{-5}\mathrm{eV}^{2}$ and for a value of $H_{0}\simeq 10^{-3}\mathrm{eV}$,
we obtain an energy density $\mathbb{T}^{\tau(MIX)}_{\tau}$ which is
compatible with the upper bound on the dark matter content of the universe. A
value of $H_{0}\simeq 10^{-3}\mathrm{eV}$ might have been reached during the
very early phases of the Universe, i.e. during the first second after the Big
Bang.

##  VI Conclusions

In this paper we have considered the quantum field theory of fermion mixing on
curved spacetime and we have computed the expectation value
$\mathbb{T}_{\mu\nu}^{(MIX)}$ of the energy-momentum tensor of fermion fields
on the flavor vacuum in a flat FLRW background. $\mathbb{T}_{\mu\nu}^{(MIX)}$
behaves as an effective energy-momentum tensor, that satisfies the Bianchi
identities (i.e. its divergence with respect to one of the two indices
vanishes) and therefore can be employed as a regular source for the Einstein
equations. In particular, it turns out that $\mathbb{T}_{\mu\nu}^{(MIX)}$ can
be interpreted as the energy-momentum tensor of a barotropic fluid.

In this picture, therefore, quantum effects of the fermion fields can be
associated, at classical level, to additional barotropic fluids whose
thermodynamical properties depend on the geometry of the spacetime other than
the specific features of the the fermion field itself. In light of the matter
field interpretation of the dark components of the Universe this result seems
quite encouraging, as it links directly quantum effects to effective classical
fluids. Thus our results show explicitly that there can be a bridge between
quantum properties of matter and dark matter/energy. However, the presence of
a connection between classical fluids and quantum effect is not automatically
sufficient alone to conclude that quantum effects are the prime cause of dark
matter/energy. In any case, our findings provide an indication that the vacuum
energy associated to neutrino mixing can represent a component of dark matter.

Here we are specifically interested in determining the form that the
energy–momentum tensor, associated to flavor vacuum, assumes in cosmological
FLRW metrics. In order to grasp the behaviour of the energy–momentum tensor of
the flavor vacuum, we have then assumed a specific form of the metric. In
doing so, we have clearly considered that the term due to the flavor vacuum is
not the dominant source of the Einstein equations, which are instead
determined by some other source that forces the specific form of the metric.

In particular, considering a de Sitter underlying metric, we have derived the
components of $\mathbb{T}_{\mu\nu}^{(MIX)}$ exactly using Bessel functions. It
turns out that at first order in the time parameter only
$\mathbb{T}_{\tau\tau}^{(MIX)}$ is different form zero. Hence, in this case
the mixing of fermion fields can be associated to a zero pressure (dust or
cold dark matter) fluid. We remark that the choice of the de Sitter metric is
dictated, primarily, by the fact that exact analytical solutions of the Dirac
equation can be found in this context. The actual identification of such a
contribution as a dark matter component requires that the study be conducted
on metrics that are adequate to the description of galaxies. Such an analysis
will be performed in a forthcoming paper, where we will also consider the
contribution due to the flavor vacuum as the dominant source in the Einstein
equations. Nevertheless, the results presented here clearly hint at a
dark–matter–like behaviour of the flavor vacuum in curved background.

##  Appendix A Orthonormality and completeness of the solution of Dirac
equation

Knowing that the determinant is $g=-C^{8}$, we can easily compute the scalar
product between the solutions (8):

$\displaystyle(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})_{\tau}$
$\displaystyle=$
$\displaystyle\int_{\Sigma_{\tau}}d^{3}xC^{3}u^{\dagger}_{\boldsymbol{p},\lambda}u_{\boldsymbol{q},\lambda^{\prime}}$
(70) $\displaystyle=$
$\displaystyle\int_{\Sigma_{\tau}}d^{3}xC^{3}e^{-i(\boldsymbol{p}-\boldsymbol{q})\cdot\boldsymbol{x}}\xi^{\dagger}_{\lambda}\xi_{\lambda^{\prime}}\left(f^{*}_{p}f_{q}+\lambda\lambda^{\prime}g^{*}_{p}g_{q}\right)$
$\displaystyle=$
$\displaystyle(2\pi)^{3}C^{3}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\delta_{\lambda,\lambda^{\prime}}\left(|f_{p}|^{2}+|g_{p}|^{2}\right)\
.$

In the last step we have used the definition of the Dirac delta function and
the orthonormality of the helicity bispinors. In a similar way it is easy to
show that
$(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=0$,
$(v_{\boldsymbol{p},\lambda},u_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=0$
and

$(v_{\boldsymbol{p},\lambda},v_{\boldsymbol{q},\lambda^{\prime}})_{\tau}=(2\pi)^{3}C^{3}\delta^{3}(\boldsymbol{p}-\boldsymbol{q})\delta_{\lambda,\lambda^{\prime}}\left(|f_{p}|^{2}+|g_{p}|^{2}\right)\
.$ (71)

Eqs.(70) and (71) suggest that we adopt the following normalization:

$|f_{p}|^{2}+|g_{p}|^{2}=\frac{1}{\left(2\pi C\right)^{3}}$ (72)

in order that the orthonormality of the solutions is satisfied. The set of
solutions (8) is complete:

$\displaystyle\sum_{\lambda}\left(u_{\boldsymbol{p},\lambda}u^{\dagger}_{\boldsymbol{p},\lambda}+v_{\boldsymbol{p},\lambda}v^{\dagger}_{\boldsymbol{p},\lambda}\right)$
$\displaystyle=$
$\displaystyle\sum_{\lambda}\left[\begin{pmatrix}f_{p}\xi_{\lambda}\\\
g_{p}\lambda\xi_{\lambda}\end{pmatrix}\begin{pmatrix}f^{*}_{p}\xi^{\dagger}_{\lambda}&g^{*}_{p}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}+\begin{pmatrix}g^{*}_{p}\xi_{\lambda}\\\
-f^{*}_{p}\lambda\xi_{\lambda}\end{pmatrix}\begin{pmatrix}g_{p}\xi^{\dagger}_{\lambda}&-f_{p}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\right]$
(73) $\displaystyle=$
$\displaystyle\sum_{\lambda}\xi_{\lambda}\xi^{\dagger}_{\lambda}\begin{pmatrix}|f_{p}|^{2}+|g_{p}|^{2}&0\\\
0&|f_{p}|^{2}+|g_{p}|^{2}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\frac{1}{\left(2\pi C\right)^{3}}\begin{pmatrix}\mathbb{I}&0\\\
0&\mathbb{I}\end{pmatrix}\ .$

In the last step we have used the completeness of the helicity basis
$\sum_{\lambda}\xi_{\lambda}\xi^{\dagger}_{\lambda}=\mathbb{I}$, with
$\mathbb{I}$ the $2\times 2$ identity matrix, and the normalization condition
(72).

##  Appendix B Properties of helicity eigenbispinors

We prove that for each $\lambda=\pm 1$ the quantity
$\xi^{\dagger}_{\lambda}(\hat{p})\sigma_{i}\xi_{\lambda}(\hat{p})$ is an odd
function of $p_{i}$, for $i=1,2,3$. In particular,
$\xi^{\dagger}_{\lambda}(\hat{p})\sigma_{i}\xi_{\lambda}(\hat{p})$ changes
sign when the momentum is reversed $\boldsymbol{p}\rightarrow-\boldsymbol{p}$.

Proof: The statement can be proven by direct calculation. For the $\lambda=+1$
case, we have

$\displaystyle\xi^{\dagger}_{+}\sigma_{1}\xi_{+}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}e^{i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})&e^{-i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\end{pmatrix}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\begin{pmatrix}e^{-i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\\\
e^{i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\end{pmatrix}$ (74)
$\displaystyle=$
$\displaystyle\cos(\frac{\theta_{p}}{2})\sin(\frac{\theta_{p}}{2})\left(e^{i\phi_{p}}+e^{-i\phi_{p}}\right)=\sin(\theta_{p})\cos(\phi_{p})=\frac{p_{1}}{p}$
$\displaystyle\xi^{\dagger}_{+}\sigma_{2}\xi_{+}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}e^{i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})&e^{-i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\end{pmatrix}\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\begin{pmatrix}e^{-i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\\\
e^{i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\end{pmatrix}$ (75)
$\displaystyle=$
$\displaystyle\cos(\frac{\theta_{p}}{2})\sin(\frac{\theta_{p}}{2})\left(ie^{-i\phi_{p}}-ie^{i\phi_{p}}\right)=\sin(\theta_{p})\sin(\phi_{p})=\frac{p_{2}}{p}$
$\displaystyle\xi^{\dagger}_{+}\sigma_{3}\xi_{+}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}e^{i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})&e^{-i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\end{pmatrix}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}e^{-i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\\\
e^{i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\end{pmatrix}$ (76)
$\displaystyle=$
$\displaystyle\cos^{2}(\frac{\theta_{p}}{2})-\sin^{2}(\frac{\theta_{p}}{2})=\cos(\theta_{p})=\frac{p_{3}}{p}$

similarly, for the $\lambda=-1$ case we have

$\displaystyle\xi^{\dagger}_{-}\sigma_{1}\xi_{-}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}e^{i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})&-e^{-i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\end{pmatrix}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\begin{pmatrix}e^{-i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\\\
-e^{i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\end{pmatrix}$ (77)
$\displaystyle=$
$\displaystyle-\cos(\frac{\theta_{p}}{2})\sin(\frac{\theta_{p}}{2})\left(e^{i\phi_{p}}+e^{-i\phi_{p}}\right)=-\sin(\theta_{p})\cos(\phi_{p})=\frac{-p_{1}}{p}$
$\displaystyle\xi^{\dagger}_{-}\sigma_{2}\xi_{-}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}e^{i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})&-e^{-i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\end{pmatrix}\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\begin{pmatrix}e^{-i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\\\
-e^{i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\end{pmatrix}$ (78)
$\displaystyle=$
$\displaystyle-\cos(\frac{\theta_{p}}{2})\sin(\frac{\theta_{p}}{2})\left(ie^{-i\phi_{p}}-ie^{i\phi_{p}}\right)=-\sin(\theta_{p})\sin(\phi_{p})=\frac{-p_{2}}{p}$
$\displaystyle\xi^{\dagger}_{-}\sigma_{3}\xi_{-}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}e^{i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})&-e^{-i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\end{pmatrix}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}e^{-i\frac{\phi_{p}}{2}}\sin(\frac{\theta_{p}}{2})\\\
-e^{i\frac{\phi_{p}}{2}}\cos(\frac{\theta_{p}}{2})\end{pmatrix}$ (79)
$\displaystyle=$
$\displaystyle-\cos^{2}(\frac{\theta_{p}}{2})+\sin^{2}(\frac{\theta_{p}}{2})=-\cos(\theta_{p})=\frac{-p_{3}}{p}\
.$

This result can be summarized in the equation

$\xi^{\dagger}_{\lambda}\boldsymbol{\sigma}\xi_{\lambda}=\frac{\lambda\boldsymbol{p}}{p}\
.$ (80)

##  Appendix C The Auxiliary Tensor

This Appendix is devoted to exploring the properties of the auxiliary tensor

${}L_{\mu\nu}(A,B)=\bar{A}\tilde{\gamma}_{\mu}(x)D_{\nu}B+\bar{A}\tilde{\gamma}_{\nu}(x)D_{\mu}B-D_{\mu}\bar{A}\tilde{\gamma}_{\nu}(x)B-D_{\nu}\bar{A}\tilde{\gamma}_{\mu}(x)B\
.$ (81)

The first property is obvious from the definition:
$L_{\mu\nu}(A,B)=L_{\nu\mu}(A,B)$ for any $A,B$. The second property can be
seen at once by comparing the definition (81) with the form of the energy-
momentum tensor (5). Since $T_{\mu\nu}$ is real, we immediately deduce that
$L_{\mu\nu}(A,A)$ is pure imaginary for each solution $A$. The third property
can be easily deduced by tracing the definition:

$L_{\mu}^{\mu}(A,B)=g^{\mu\nu}L_{\mu\nu}(A,B)=2\left(\bar{A}\tilde{\gamma}^{\mu}(x)D_{\mu}B-D_{\mu}\bar{A}\tilde{\gamma}_{\mu}(x)B\right)=-4im\bar{A}B$
(82)

where we have used the Dirac equation and its adjoint in the last step. Next
we give a more explicit form for $L_{\tau\tau}(A,B)$. Given the expression of
the $\tilde{\gamma}_{\mu}$ matrices we have

$\displaystyle L_{\tau\tau}(A,B)$ $\displaystyle=$ $\displaystyle
2\left(\bar{A}\tilde{\gamma}_{\tau}(x)D_{\tau}B-D_{\tau}\bar{A}\tilde{\gamma}_{\tau}(x)B\right)=2C\left(\bar{A}\gamma^{0}\partial_{\tau}B-\partial_{\tau}\bar{A}\gamma^{0}B\right)$
(83) $\displaystyle=$ $\displaystyle
2C\left(A^{\dagger}\partial_{\tau}B-\partial_{\tau}A^{\dagger}B\right)\ .$

In particular we are interested in the auxiliary tensors
$L_{\tau\tau}(u(v),u(v))$ and the traces $L^{\mu}_{\mu}(u(v),u(v))$computed on
the solutions (8). For the traces we have

$\displaystyle
L^{\mu}_{\mu}(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})$
$\displaystyle=$
$\displaystyle-4im\bar{u}_{\boldsymbol{p},\lambda}u_{\boldsymbol{p},\lambda}=-4imu^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}u_{\boldsymbol{p},\lambda}$
(84) $\displaystyle=$
$\displaystyle-4im\begin{pmatrix}f^{*}_{p}\xi^{\dagger}_{\lambda}&g^{*}_{p}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}\mathbb{I}&0\\\
0&-\mathbb{I}\end{pmatrix}\begin{pmatrix}f_{p}\xi_{\lambda}\\\
g_{p}\lambda\xi_{\lambda}\end{pmatrix}$ $\displaystyle=$
$\displaystyle-4im(|f_{p}|^{2}-|g_{p}|^{2})=-4imC^{-3}\left(|\phi_{p}|^{2}-|\gamma_{p}|^{2}\right)$
$\displaystyle
L^{\mu}_{\mu}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})$
$\displaystyle=$
$\displaystyle-4im\bar{u}_{\boldsymbol{p},\lambda}v_{\boldsymbol{p},\lambda}=-4imu^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}v_{\boldsymbol{p},\lambda}$
(85) $\displaystyle=$
$\displaystyle-4im\begin{pmatrix}f^{*}_{p}\xi^{\dagger}_{\lambda}&g^{*}_{p}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}\mathbb{I}&0\\\
0&-\mathbb{I}\end{pmatrix}\begin{pmatrix}g_{p}^{*}\xi_{\lambda}\\\
-f^{*}_{p}\lambda\xi_{\lambda}\end{pmatrix}$ $\displaystyle=$
$\displaystyle-8imf^{*}_{p}g^{*}_{p}=-8imC^{-3}\phi^{*}_{p}\gamma^{*}_{p}$
$\displaystyle
L^{\mu}_{\mu}(v_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})$
$\displaystyle=$
$\displaystyle-4im\bar{v}_{\boldsymbol{p},\lambda}u_{\boldsymbol{p},\lambda}=-4imv^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}u_{\boldsymbol{p},\lambda}$
(86) $\displaystyle=$
$\displaystyle-4im\begin{pmatrix}g_{p}\xi^{\dagger}_{\lambda}&-f_{p}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}\mathbb{I}&0\\\
0&-\mathbb{I}\end{pmatrix}\begin{pmatrix}f_{p}\xi_{\lambda}\\\
g_{p}\lambda\xi_{\lambda}\end{pmatrix}$ $\displaystyle=$
$\displaystyle-8imf_{p}g_{p}=-8imC^{-3}\phi_{p}\gamma_{p}=-\left(L^{\mu}_{\mu}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})\right)^{*}$
$\displaystyle
L^{\mu}_{\mu}(v_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})$
$\displaystyle=$
$\displaystyle-4im\bar{v}_{\boldsymbol{p},\lambda}v_{\boldsymbol{p},\lambda}=-4imv^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}v_{\boldsymbol{p},\lambda}$
(87) $\displaystyle=$
$\displaystyle-4im\begin{pmatrix}g_{p}\xi^{\dagger}_{\lambda}&-f_{p}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}\mathbb{I}&0\\\
0&-\mathbb{I}\end{pmatrix}\begin{pmatrix}g^{*}_{p}\xi_{\lambda}\\\
-f^{*}_{p}\lambda\xi_{\lambda}\end{pmatrix}$ $\displaystyle=$
$\displaystyle-4im(|g_{p}|^{2}-|f_{p}|^{2})=-4imC^{-3}\left(|\gamma_{p}|^{2}-|\phi_{p}|^{2}\right)$
$\displaystyle=$
$\displaystyle-L^{\mu}_{\mu}(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})\
.$

We then compute the $\tau\tau$ components

$\displaystyle
L_{\tau\tau}(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})$
$\displaystyle=$ $\displaystyle
2C\left[u^{\dagger}_{\boldsymbol{p},\lambda}\partial_{\tau}u_{\boldsymbol{p},\lambda}-\partial_{\tau}u^{\dagger}_{\boldsymbol{p},\lambda}u_{\boldsymbol{p},\lambda}\right]$
(88) $\displaystyle=$ $\displaystyle
2C\left[\begin{pmatrix}f_{p}^{*}\xi^{\dagger}_{\lambda}&g_{p}^{*}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}\partial_{\tau}f_{p}\xi_{\lambda}\\\
\partial_{\tau}g_{p}\lambda\xi_{\lambda}\end{pmatrix}-\begin{pmatrix}\partial_{\tau}f_{p}^{*}\xi^{\dagger}_{\lambda}&\partial_{\tau}g_{p}^{*}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}f_{p}\xi_{\lambda}\\\
g_{p}\lambda\xi_{\lambda}\end{pmatrix}\right]$ $\displaystyle=$ $\displaystyle
2C\left[f^{*}_{p}\partial_{\tau}f_{p}-\partial_{\tau}f^{*}_{p}f_{p}+g^{*}_{p}\partial_{\tau}g_{p}-\partial_{\tau}g^{*}_{p}g_{p}\right]$
$\displaystyle=$ $\displaystyle
2C^{-2}\left[\phi^{*}_{p}\partial_{\tau}\phi_{p}-\partial_{\tau}\phi^{*}_{p}\phi_{p}+\gamma^{*}_{p}\partial_{\tau}\gamma_{p}-\partial_{\tau}\gamma^{*}_{p}\gamma_{p}\right]$
$\displaystyle
L_{\tau\tau}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})$
$\displaystyle=$ $\displaystyle
2C\left[u^{\dagger}_{\boldsymbol{p},\lambda}\partial_{\tau}v_{\boldsymbol{p},\lambda}-\partial_{\tau}u^{\dagger}_{\boldsymbol{p},\lambda}v_{\boldsymbol{p},\lambda}\right]$
(89) $\displaystyle=$ $\displaystyle
2C\left[\begin{pmatrix}f_{p}^{*}\xi^{\dagger}_{\lambda}&g_{p}^{*}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}\partial_{\tau}g_{p}^{*}\xi_{\lambda}\\\
-\partial_{\tau}f_{p}^{*}\lambda\xi_{\lambda}\end{pmatrix}-\begin{pmatrix}\partial_{\tau}f_{p}^{*}\xi^{\dagger}_{\lambda}&\partial_{\tau}g_{p}^{*}\lambda\xi^{\dagger}_{\lambda}\end{pmatrix}\begin{pmatrix}g_{p}^{*}\xi_{\lambda}\\\
-f^{*}_{p}\lambda\xi_{\lambda}\end{pmatrix}\right]$ $\displaystyle=$
$\displaystyle
2C\left[f^{*}_{p}\partial_{\tau}g_{p}^{*}-\partial_{\tau}f^{*}_{p}g_{p}^{*}-g^{*}_{p}\partial_{\tau}f_{p}^{*}+\partial_{\tau}g^{*}_{p}f^{*}_{p}\right]$
$\displaystyle=$ $\displaystyle
4C^{-2}\left[\phi^{*}_{p}\partial_{\tau}\gamma^{*}_{p}-\gamma^{*}_{p}\partial_{\tau}\phi_{p}^{*}\right]$
$\displaystyle
L_{\tau\tau}(v_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})$
$\displaystyle=$
$\displaystyle-\left(L_{\tau\tau}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})\right)^{*}$
(90) $\displaystyle
L_{\tau\tau}(v_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})$
$\displaystyle=$ $\displaystyle-
L_{\tau\tau}(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})\ .$ (91)

Notice that as a consequence
$L_{\tau\tau}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})=0$ if and
only if $\gamma_{p}\propto\phi_{p}$. This is the case, for instance, in
Minkowski spacetime, where $\phi_{p}\propto\gamma_{p}\propto
e^{-i\omega_{p}t}$. For a general $C$, the proportionality does not hold, and
this has dramatic consequences on the vacuum of the theory (eventually leading
to particle creation).

In addition one can show that $L_{\tau i}(u(v),u(v))$ is an odd function of
the momentum $\boldsymbol{p}$ and that $L_{ij}(u(v),u(v))$ is an odd function
with respect to both $p_{i}$ and $p_{j}$. The proof requires a lengthy but
straightforward calculation:

$\displaystyle L_{\tau
i}(u_{\boldsymbol{p},\lambda},u_{\boldsymbol{p},\lambda})$ $\displaystyle=$
$\displaystyle\bar{u}_{\boldsymbol{p},\lambda}\tilde{\gamma}_{\tau}D_{i}u_{\boldsymbol{p},\lambda}+\bar{u}_{\boldsymbol{p},\lambda}\tilde{\gamma}_{i}D_{\tau}u_{\boldsymbol{p},\lambda}-D_{i}\bar{u}_{\boldsymbol{p},\lambda}\tilde{\gamma}_{\tau}u_{\boldsymbol{p},\lambda}-D_{\tau}\bar{u}_{\boldsymbol{p},\lambda}\tilde{\gamma}_{i}u_{\boldsymbol{p},\lambda}$
(92) $\displaystyle=$ $\displaystyle
Cu^{\dagger}_{\boldsymbol{p},\lambda}\left(\partial_{i}+\frac{1}{8}\omega_{i}^{A,B}[\gamma_{A},\gamma_{B}]\right)u_{\boldsymbol{p},\lambda}-Cu^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}\gamma^{i}\partial_{\tau}u_{\boldsymbol{p},\lambda}$
$\displaystyle-$ $\displaystyle
C\left(\partial_{i}u^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}-u^{\dagger}_{\boldsymbol{p},\lambda}\frac{\gamma^{0}}{8}\omega_{i}^{A,B}[\gamma_{A},\gamma_{B}]\right)\gamma^{0}u_{\boldsymbol{p},\lambda}+C\partial_{\tau}u^{\dagger}_{\boldsymbol{p},\lambda}\gamma^{0}\gamma^{i}u_{\boldsymbol{p},\lambda}$
$\displaystyle=$ $\displaystyle
Cu^{\dagger}_{\boldsymbol{p},\lambda}\left(ip_{i}+\frac{\partial_{\tau}C}{2C}\begin{pmatrix}0&\sigma_{i}\\\
\sigma_{i}&0\end{pmatrix}\right)u_{\boldsymbol{p},\lambda}-Cu^{\dagger}_{\boldsymbol{p},\lambda}\begin{pmatrix}0&\sigma_{i}\\\
\sigma_{i}&0\end{pmatrix}\partial_{\tau}u_{\boldsymbol{p},\lambda}$
$\displaystyle-$ $\displaystyle
Cu^{\dagger}_{\boldsymbol{p},\lambda}\left(-ip_{i}+\frac{\partial_{\tau}C}{2C}\begin{pmatrix}0&\sigma_{i}\\\
\sigma_{i}&0\end{pmatrix}\right)u_{\boldsymbol{p},\lambda}+C\partial_{\tau}u^{\dagger}_{\boldsymbol{p},\lambda}\begin{pmatrix}0&\sigma_{i}\\\
\sigma_{i}&0\end{pmatrix}u_{\boldsymbol{p},\lambda}$ $\displaystyle=$
$\displaystyle
2ip_{i}Cu^{\dagger}_{\boldsymbol{p},\lambda}u_{\boldsymbol{p},\lambda}+C\left[\partial_{\tau}u^{\dagger}_{\boldsymbol{p},\lambda}\begin{pmatrix}0&\sigma_{i}\\\
\sigma_{i}&0\end{pmatrix}u_{\boldsymbol{p},\lambda}-u^{\dagger}_{\boldsymbol{p},\lambda}\begin{pmatrix}0&\sigma_{i}\\\
\sigma_{i}&0\end{pmatrix}\partial_{\tau}u_{\boldsymbol{p},\lambda}\right]$
$\displaystyle=$ $\displaystyle
2ip_{i}C\left[|f_{p}|^{2}+|g_{p}|^{2}\right]+C\lambda\left[\partial_{\tau}f_{p}^{*}g_{p}+\partial_{\tau}g_{p}^{*}f_{p}-f_{p}^{*}\partial_{\tau}g_{p}-g_{p}^{*}\partial_{\tau}f_{p}\right]\xi^{\dagger}_{\lambda}\sigma_{i}\xi_{\lambda}$
$\displaystyle=$ $\displaystyle
p_{i}\left\\{\frac{i}{\pi^{3}C^{2}}+\frac{\lambda}{p}\left[C\left(\partial_{\tau}f_{p}^{*}g_{p}+\partial_{\tau}g_{p}^{*}f_{p}-c.c.\right)\right]\right\\}\
.$

In the last step, we have made use of the normalization condition (72) and the
property

(80). As it is evident from Eq.(92), each $\tilde{\gamma}_{i}$ factor and each
spatial derivative $\partial_{i}$ brings along a factor $p_{i}$. Then for each
$a,b=u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda}$ one has

$L_{\tau i}(a,b)=p_{i}h_{a,b}(p)\ ,$ (93)

with $h_{a,b}(p)$ a function of the modulus $p$ alone. In particular

$L_{\tau i}(u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda})=0\ .$ (94)

Similarly, for each
$a,b=u_{\boldsymbol{p},\lambda},v_{\boldsymbol{p},\lambda}$

$L_{ij}(a,b)=p_{i}p_{j}l_{a,b}(p)\ ,$ (95)

with $l_{a,b}(p)$ a function of the modulus $p$ alone.

##  Appendix D Canonical anticommutation relations

Demonstration of the canonical anticommutation relations (20)

$\displaystyle\\{\psi_{A}(\tau,\boldsymbol{x}),\pi_{\psi
B}(\tau,\boldsymbol{x^{\prime}})\\}$ $\displaystyle=$ $\displaystyle
iC^{3}\sum_{\lambda,\lambda^{\prime}}\int d^{3}p\int
d^{3}q\left[\\{A_{\boldsymbol{p},\lambda},A^{\dagger}_{\boldsymbol{q},\lambda^{\prime}}\\}\left(u_{\boldsymbol{p},\lambda}u^{\dagger}_{\boldsymbol{q},\lambda^{\prime}}\right)_{AB}+\\{B^{\dagger}_{-\boldsymbol{p},\lambda},B_{-\boldsymbol{q},\lambda^{\prime}}\\}\left(v_{\boldsymbol{p},\lambda}v^{\dagger}_{\boldsymbol{q},\lambda^{\prime}}\right)_{AB}\right]$
(96) $\displaystyle=$ $\displaystyle iC^{3}\sum_{\lambda}\int
d^{3}p\left[\left(u_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{x})u^{\dagger}_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{x^{\prime}})\right)_{AB}+\left(v_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{x})v^{\dagger}_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{x^{\prime}})\right)_{AB}\right]$
$\displaystyle=$ $\displaystyle iC^{3}\sum_{\lambda}\int
d^{3}pe^{i\boldsymbol{p}\cdot(\boldsymbol{x}-\boldsymbol{x^{\prime}})}\left[\left(u_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{0})u^{\dagger}_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{0})\right)_{AB}+\left(v_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{0})v^{\dagger}_{\boldsymbol{p},\lambda}(\tau,\boldsymbol{0})\right)_{AB}\right]$
$\displaystyle=$ $\displaystyle iC^{3}\frac{1}{(2\pi)^{3}C^{3}}\delta_{AB}\int
d^{3}pe^{i\boldsymbol{p}\cdot(\boldsymbol{x}-\boldsymbol{x^{\prime}})}$
$\displaystyle=$ $\displaystyle
i\delta_{AB}\delta^{3}(\boldsymbol{x}-\boldsymbol{x^{\prime}})\ .$

In the fourth step we have made use of the completeness relation (Eq.(73)),
writing the $4\times 4$ identity matrix explicitly as $\delta_{AB}$.

##  Appendix E Conservation of the energy–momentum tensor

In this appendix we prove explicitly the covariant conservation of the energy
momentum tensor associated to the flavor vacuum. We show that

$\nabla_{\mu}\mathbb{T}^{\mu\nu}=0$ (97)

with $\nabla_{\mu}$ denoting the covariant derivative. There is no need here
to distinguish between $\mathbb{T}_{\mu\nu}^{(MIX)}$ and
$\mathbb{T}_{\mu\nu}^{(N)}$, since both satisfy eq. (97), and so does the full
energy-momentum tensor. Preliminarly we derive the connection coefficients for
the metric of eq. (1). It is easy to see from the definition that the only non
vanishing coefficients are

$\Gamma^{\tau}_{\tau\tau}=\Gamma^{\tau}_{ii}=\Gamma^{i}_{\tau
i}=\Gamma^{i}_{i\tau}=\frac{\dot{C}}{C}\ .$ (98)

Here the dot denotes the derivative with respect to conformal time $\tau$ and
no sum is intended over repeated indices. Notice that the coefficients depend
only on $\tau$. In terms of the connection coefficients, the covariant
divergence reads

$\nabla_{\mu}\mathbb{T}^{\mu\nu}=\partial_{\mu}\mathbb{T}^{\mu\nu}+\Gamma^{\mu}_{\mu\sigma}\mathbb{T}^{\sigma\nu}+\Gamma^{\nu}_{\mu\sigma}\mathbb{T}^{\mu\sigma}\
.$ (99)

  * •

($\nu=i$) For $\nu=i$, with $i=1,2,3$ equation (99) becomes

$\nabla_{\mu}\mathbb{T}^{\mu i}=\partial_{\mu}\mathbb{T}^{\mu
i}+\Gamma^{\mu}_{\mu\sigma}\mathbb{T}^{\sigma
i}+\Gamma^{i}_{\mu\sigma}\mathbb{T}^{\mu\sigma}\ .$ (100)

From the diagonality of $\mathbb{T}^{\mu\nu}$ proved above, we can write

$\nabla_{\mu}\mathbb{T}^{\mu
i}=\partial_{i}\mathbb{T}^{ii}+\sum_{\mu}\Gamma^{\mu}_{\mu
i}\mathbb{T}^{ii}+\sum_{\mu}\Gamma^{i}_{\mu\mu}\mathbb{T}^{\mu\mu}\ ,$ (101)

where no sum is intended over repeated indices and the summations are written
out explicitly to avoid confusion. The first term on the right hand side of
eq. (101) is zero, since $\mathbb{T}^{\mu\nu}$ depends only on $\tau$.
Similarly, from eq. (98) we know that $\Gamma^{\mu}_{\mu
i}=0=\Gamma^{i}_{\mu\mu}$ for each $\mu=0,1,2,3$ and each $i=1,2,3$, so that
also the second and the third term on the right hand side of eq. (101) vanish.
Then

$\nabla_{\mu}\mathbb{T}^{\mu i}=0\ \ \ \ \forall i\ .$ (102)

  * •

($\nu=\tau$) Only a slightly longer calculation is needed to prove the
statement for $\nu=\tau$. Starting from equation (99) we have

$\displaystyle\nabla_{\mu}\mathbb{T}^{\mu\tau}$ $\displaystyle=$
$\displaystyle\partial_{\mu}\mathbb{T}^{\mu\tau}+\Gamma^{\mu}_{\mu\sigma}\Gamma^{\sigma\tau}+\Gamma^{\tau}_{\mu\sigma}\mathbb{T}^{\mu\sigma}$
(103) $\displaystyle=$
$\displaystyle\partial_{\tau}\mathbb{T}^{\tau\tau}+\left(\Gamma^{\tau}_{\tau\tau}+\sum_{i}\Gamma^{i}_{i\tau}\right)\mathbb{T}^{\tau\tau}+\Gamma^{\tau}_{\tau\tau}\mathbb{T}^{\tau\tau}+\sum_{i}\Gamma^{\tau}_{ii}\mathbb{T}^{ii}$
$\displaystyle=$
$\displaystyle\partial_{\tau}\mathbb{T}^{\tau\tau}+5\Gamma^{\tau}_{\tau\tau}\mathbb{T}^{\tau\tau}+3\Gamma^{\tau}_{\tau\tau}\mathbb{T}^{ii}\
,$

where we have used the diagonality of $\mathbb{T}^{\mu\nu}$ and eqs. (98). For
our purposes it is convenient to rewrite eq. (103) in terms of
$\mathbb{T}_{\tau\tau}$ and the trace $\mathbb{T}^{\mu}_{\mu}$. To this end we
employ eqs. (47), (98) and lower the indices through the metric of eq. (1),
obtaining:

$\nabla_{\mu}\mathbb{T}^{\mu\tau}=\partial_{\tau}\left(C^{-4}\mathbb{T}_{\tau\tau}\right)+6C^{-5}\dot{C}\mathbb{T}_{\tau\tau}-C^{-3}\dot{C}\mathbb{T}^{\mu}_{\mu}\
.$ (104)

From equation (33), we know that each of the terms above is the integral of
the auxiliary tensor components $L_{\tau\tau},L^{\mu}_{\mu}$ weighted by
$\tau$-independent coefficients (recall that the Bogoliubov coefficients are
evaluated at a fixed time $\tau_{0}$). It is therefore sufficient to prove
that

$\partial_{\tau}\left(C^{-4}L_{\tau\tau}(a,b)\right)+6C^{-5}\dot{C}L_{\tau\tau}(a,b)-C^{-3}\dot{C}L^{\mu}_{\mu}(a,b)=0$
(105)

for each $a,b=u_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j}$, to
show that the divergence (104) vanishes. To this end we need the second order
equations

$\displaystyle\partial^{2}_{\tau}\phi_{p;j}$ $\displaystyle=$
$\displaystyle-\left(im_{j}\dot{C}+p^{2}+m_{j}^{2}C^{2}\right)\phi_{p;j}$
$\displaystyle\partial^{2}_{\tau}\gamma_{p;j}$ $\displaystyle=$
$\displaystyle-\left(-im_{j}\dot{C}+p^{2}+m_{j}^{2}C^{2}\right)\gamma_{p;j}\
.$ (106)

The first of this equations is simply eq. (17), and the second can be likewise
deduced from the system (15). For
$L_{\mu\nu}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})$ we
have, using the properties of the auxiliary tensor

$\displaystyle\partial_{\tau}\left(C^{-4}L_{\tau\tau}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})\right)+6C^{-5}\dot{C}L_{\tau\tau}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})-C^{-3}\dot{C}L_{\mu}^{\mu}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})=$
$\displaystyle
2\partial_{\tau}\left[C^{-6}\left(\phi_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\phi_{p;j}+\gamma_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\gamma_{p;j}\right)\right]+12C^{-7}\dot{C}\left(\phi_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\phi_{p;j}+\gamma_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\gamma_{p;j}\right)$
$\displaystyle+4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)=$
$\displaystyle
2C^{-6}\partial_{\tau}\left(\phi_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\phi_{p;j}+\gamma_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\gamma_{p;j}\right)+4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)=$
$\displaystyle
2C^{-6}\left(\phi_{p;j}^{*}\partial^{2}_{\tau}\phi_{p;j}-\phi_{p;j}\partial_{\tau}^{2}\phi_{p;j}^{*}+\gamma_{p;j}^{*}\partial^{2}_{\tau}\gamma_{p;j}-\gamma_{p;j}\partial_{\tau}^{2}\gamma_{p;j}^{*}\right)+4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)=$
$\displaystyle
2C^{-6}\Big{\\{}\phi_{p;j}^{*}\left[-\left(im_{j}\dot{C}+p^{2}+m^{2}C^{2}\right)\right]\phi_{p;j}-\phi_{p;j}^{*}\left[-\left(-im_{j}\dot{C}+p^{2}+m^{2}C^{2}\right)\right]\phi_{p;j}$
$\displaystyle+\gamma_{p;j}^{*}\left[-\left(-im_{j}\dot{C}+p^{2}+m^{2}C^{2}\right)\right]\gamma_{p;j}-\phi_{p;j}^{*}\left[-\left(im_{j}\dot{C}+p^{2}+m^{2}C^{2}\right)\right]\gamma_{p;j}\Big{\\}}$
$\displaystyle+4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)=$
$\displaystyle
2C^{-6}\left\\{|\phi_{p;j}|^{2}\left(-2im_{j}\dot{C}\right)+|\gamma_{p;j}|^{2}\left(2im_{j}\dot{C}\right)\right\\}4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)=$
$\displaystyle-4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)+4im_{j}C^{-6}\dot{C}\left(|\phi_{p;j}|^{2}-|\gamma_{p;j}|^{2}\right)=0\
.$

In the fourth step we have used the equations (• ‣ E) and their complex
conjugates, and the $\tau$ argument of the functions has been suppressed for
notational simplicity. Note that from the properties
$L_{\tau\tau}(v_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})=-L_{\tau\tau}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})$
and
$L_{\mu}^{\mu}(v_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})=-L_{\mu}^{\mu}(u_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})$
the same relation is satisfied by the components of
$L_{\mu\nu}(v_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})$.
Similarly one has

$\displaystyle\partial_{\tau}\left(C^{-4}L_{\tau\tau}(u_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})\right)+6C^{-5}\dot{C}L_{\tau\tau}(u_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})-C^{-3}\dot{C}L_{\mu}^{\mu}(u_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})=$
$\displaystyle\partial_{\tau}\left(4C^{-6}\phi_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\gamma_{p;j}^{*}\right)+24C^{-7}\dot{C}\left(4C^{-6}\phi_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\gamma_{p;j}^{*}\right)+8im_{j}C^{-6}\dot{C}\phi_{p;j}^{*}\gamma_{p;j}^{*}=$
$\displaystyle
4C^{-6}\partial_{\tau}\left(C^{-6}\phi_{p;j}^{*}\overleftrightarrow{\partial_{\tau}}\gamma_{p;j}^{*}\right)+8im_{j}C^{-6}\dot{C}\phi_{p;j}^{*}\gamma_{p;j}^{*}=$
$\displaystyle
4\left[\phi_{p;j}^{*}\partial^{2}_{\tau}\gamma^{*}_{p;j}-\gamma^{*}_{p;j}\partial^{2}_{\tau}\phi_{p;j}^{*}\right]+8im_{j}C^{-6}\dot{C}\phi_{p;j}^{*}\gamma_{p;j}^{*}=$
$\displaystyle
4C^{-6}\left\\{\phi_{p;j}^{*}\left[-\left(im_{j}\dot{C}+p^{2}+m^{2}C^{2}\right)\right]\gamma_{p;j}^{*}-\phi_{p;j}^{*}\left[-\left(-im_{j}\dot{C}+p^{2}+m^{2}C^{2}\right)\right]\gamma_{p;j}^{*}\right\\}$
$\displaystyle+8im_{j}C^{-6}\dot{C}\phi_{p;j}^{*}\gamma_{p;j}^{*}=$
$\displaystyle-8im_{j}C^{-6}\dot{C}\phi_{p;j}^{*}\gamma_{p;j}^{*}++8im_{j}C^{-6}\dot{C}\phi_{p;j}^{*}\gamma_{p;j}^{*}=0\
.$

In the fourth step we have used the complex conjugates of eqs. (• ‣ E).
Finally, due to the properties
$L_{\tau\tau}(v_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})=-L^{*}_{\tau\tau}(u_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})$
and
$L_{\mu}^{\mu}(v_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})=-L_{\mu}^{\mu*}(u_{\boldsymbol{p},\lambda;j},v_{\boldsymbol{p},\lambda;j})$,
the same relation is satisfied by the components of
$L_{\mu\nu}(v_{\boldsymbol{p},\lambda;j},u_{\boldsymbol{p},\lambda;j})$. This
is sufficient to prove the statement for $\nu=\tau$.

## References

  * [1] P. de Bernardis et al., Nature, 404, 955 (2000).
  * [2] D.N. Spergel et al. ApJS, 148, 175 (2003).
  * [3] S. Dodelson et al., ApJ, 572, 140 (2002).
  * [4] A.S. Szalay et al., ApJ, 591, 1 (2003).
  * [5] Durrer, R., Maartens, R. (2008). General Relativity and Gravitation, 40(2-3), 301-328.
  * [6] P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003).
  * [7] A.G. Riess et al., ApJ, 607, 665 (2004).
  * [8] V. Trimble, Annual Review of Astronomy and Astrophysics 25.1, pp. 425–472 (1987).
  * [9] E. Corbelli and P. Salucci, Monthly Notices of the Royal Astronomical Society 311.2, pp. 441-447 (Jan. 2000).
  * [10] V. C. Rubin, Jr. Ford W. K., Astrophysical Journal 159, p. 379 (Feb. 1970).
  * [11] V. C. Rubin, Jr. Ford W. K. and N. Thonnard, Astrophysical Journal 238, pp. 471–487 (June 1980).
  * [12] D. Clowe et al., ApJ 649, 2 (2006).
  * [13] N. Aghanim et al., A & A 641, A6 (2020).
  * [14] A. Joyce et al., Phys. Rep. 568, pp. 1–98 (2015).
  * [15] G. Bertone, D. Hooper, J. Silk, Phys. Rep. 405, issues 5–6, pp. 279–390 (2005).
  * [16] P. F. de Salas and A. Widmark, Rep. Prog. Phys. 84, 104901 (2021).
  * [17] G. Bertone and D. Hooper, “History of dark matter”, Rev. Mod. Phys. 90, 045002 (2018).
  * [18] L. Bergström, Rep. Prog. Phys. 63, 793 (2000).
  * [19] L. Roszkowski, E. M. Sessolo and S. Trojanowski, Rep. Prog. Phys. 81, 066201 (2018).
  * [20] M. Klasen, M. Pohl and G. Sigl, Prog. Part. Nucl. Phys. 85, pp. 1–32 (2015).
  * [21] J. Conrad and O. Reimer, Nat. Phys. 13, 224 (2017).
  * [22] M. Schumann, J. Phys. G: Nucl. Part. Phys. 46, 103003 (2019).
  * [23] A. Berlin, Phys. Rev. Lett. 117, 231801 (2016).
  * [24] L. Necib, M. Lisanti and V. Belokurov, Astrophys. J 874, 3 (2019).
  * [25] F. Mayet et al., Phys. Rep. 627, 1 (2016).
  * [26] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, pp. 1791-1797 (1977).
  * [27] A. Arvanitaki et al., Phys. Rev. D 81 (2010)
  * [28] G. Raffelt, Stars as Laboratories for Fundamental Physics, University of Chicago Press, Chicago, ill, USA, 1996.
  * [29] G. G. Raffelt, “Astrophysical axion bounds”, in Axions, Lecture Notes in Physics Vol. 741, pp. 51–71, Springer, Berlin, Germany, 2008.
  * [30] A. Capolupo, G. Lambiase, A. Quaranta, S. M. Giampaolo, Phys. Lett. B 804, 135407 (2020).
  * [31] I. De Martino, T. Broadhurst, S. H. H. Tye, T. Chiueh, H.Y. Schive, R. Lazkoz, Phys. Rev. Lett. 119 (22), 221103 (2017).
  * [32] A. Capolupo, I. De Martino, G. Lambiase and An. Stabile, Phys. Lett. B 790, pp. 427–435 (2019).
  * [33] I. De Martino, T. Broadhurst, S. H. H. Tye, T. Chiueh, H. Y. Schive, R. Lazkoz, Galaxies 6 (1), 10 (2018).
  * [34] A. Capolupo, S. M. Giampaolo, A. Quaranta, “Neutron interferometry and axionlike particles”, arXiv:2102.13206v1 [hep-ph] (2021).
  * [35] S. Aune et al., Phys. Rev. Lett. 107, no. 26, article ID 261302 (2011).
  * [36] A. Capolupo, G. Lambiase, G. Vitiello, Adv. High En. Phys., 2015, 826051, 7 pages (2015).
  * [37] G. Jungman, M. Kamionkowski and K. Griest, Phys. Rep. 267.5, pp. 195–373 (1996).
  * [38] S. Clesse and J. Garcia–Bellido, Phys. Dark Univ. 22, pp. 137–146 (2018).
  * [39] A. Capolupo, Adv. High En. Phys. 2018, 9840351 (2018).
  * [40] A. Capolupo, Adv. High En. Phys. 2016, 8089142 (2016).
  * [41] A. Capolupo, S. Capozziello, G. Vitiello, Phys. Lett. A 373, pp. 601–610 (2009).
  * [42] A. Capolupo, S. Capozziello, G. Vitiello, Phys. Lett. A 363, 53 (2007).
  * [43] M. Blasone, A. Capolupo, S. Capozziello, S. Carloni, G. Vitiello, Phys. Lett. A 323, pp. 182–189 (2004).
  * [44] W. Buchmüller, Neutrinos, Grand Unification and Leptogenesis, arXiv:hep-ph/0204288v2 (2002).
  * [45] M. Tanabashi et al. (Particle Data Group), Neutrinos in Cosmology, Phys. Rev. D 98, 030001 (2018).
  * [46] M.G. Betti et al. (PTOLEMY collaboration), JCAP 07, 047 (2019).
  * [47] D. B. Kaplan, A. E. Nelson, N. Weiner, Phys. Rev. Lett. 93, 091801 (2004); R. Fardon, A. E. Nelson, N. Weiner, JCAP 10, 005 (2004).
  * [48] A. Capolupo, G. Lambiase, A. Quaranta Phys. Rev. D, 101, 095022 (2020).
  * [49] Y. Grossman and H. J. Lipkin, Phys. Rev. D 55, 2760 (1997).
  * [50] D. Piriz, M. Roy and J. Wudka, Phys. Rev. D 54, 1587 (1996).
  * [51] C. Y. Cardall and G. M. Fuller, Phys. Rev. D 55, 7960 (1997).
  * [52] L. Buoninfante, G. G. Luciano, L. Petruzziello and L. Smaldone, Phys. Rev. D 101, 024016 (2020).
  * [53] L. Parker, Phys. Rev. 183, 1057 (1969) and Phys. Rev. D 3, 346 (1971).
  * [54] N. D. Birrell and P. C. Davies, “Quantum fields in curved space”, Cambridge University Press (1982).
  * [55] L. H. Ford, “Quantum field theory in curved space–time”, 9th Jorge Andre Swieca Summer School: Particles and Fields, arXiv gr-qc/9707062, pp. 345–388 (1997)
  * [56] L. Parker and D. Toms, “Quantum field theory in curved spacetime: quantized fields and gravity”, Cambridge University Press (2009).
  * [57] M. Abramowitz and I. A. Stegun, Handbook Of Mathematical Functions, Dover Publications Inc., New York, ISBN-13: 978-0-486-61272-0.