# Distinguishing Dark Energy Models with Neutrino Oscillations

Ali Rida Khalifeh<EMAIL_ADDRESS>Raul Jimenez<EMAIL_ADDRESS>ICC,
University of Barcelona, Marti i Franques, 1, E-08028 Barcelona, Spain.  Dept.
de Fisica Cuantica y Astrofisica, University of Barcelona, Marti i Franques 1,
E-08028 Barcelona, Spain.  ICREA, Pg. Lluis Companys 23, Barcelona, E-08010,
Spain.

###### Abstract

Dark Energy models are numerous and distinguishing between them is becoming
difficult. However, using distinct observational probes can ease this quest
and gives better assessment to the nature of Dark energy. To this end, the
plausibility of neutrino oscillations to be a probe of Dark Energy models is
investigated. First, a generalized formalism of neutrino (spinor field)
interaction with a classical scalar field in curved space-time is presented.
This formalism is then applied to two classes of Dark Energy models in a flat
Friedman-Lemaître-Robertson-Walker metric: a Cosmological Constant and scalar
field Dark Energy coupled to neutrinos. By looking at the neutrino oscillation
probability’s evolution with redshift, these models can be distinguished, for
certain neutrino and scalar field coupling properties. This evolution could be
traced by neutrino flux measurements in future underground, terrestrial or
extraterrestrial, neutrino telescopes which would assess probing Dark Energy
models with this technique.

##  1 Introduction

Since the discovery of the accelerated expansion of the Universe [1, 2], one
of the most interesting open questions in Astrophysics and Cosmology is to
understand if Dark Energy(DE) is dynamic, or instead strictly a constant.
Indeed, if DE was shown to be dynamical, this would be a major revolution, as
it would indicate a great deal of new physics. However, recent observational
constraints indicate that DE is consistent with a cosmological constant, with
a few percent uncertainty [3]. Current and upcoming cosmological surveys, such
as DESI [4] and Euclid [5] will decrease this level of uncertainty to the %
level [6]. Nevertheless, theoretical arguments have been presented over the
years in favor of a dynamical DE, due to fundamental issues accompanied by a
constant one, such as the coincidence problem [7], see also [8, 9, 10].

In order to lift this dilemma, one could combine several probes and techniques
to constraint DE models. In addition to the already mentioned probes, as well
as Gravitational Waves surveys [11, 12, 13], looking at neutrinos could open a
new window to the nature of DE. Several Cosmological probes have been used to
constrain neutrino properties in the context of a flat Friedman-Lemaître-
Robertson-Walker(FLRW) Universe [14, 15, 16, 17, 18, 19]. However, neutrinos
have been mostly considered as classical point particles, rather than quantum
spinor fields traveling in curved spacetime, which could provide novel
insights for both neutrinos and DE. An interaction between spinor and scalar
fields is sensible to our spacetime’s curved geometry, which could leave
observational imprints. This shows the advantage of using Cosmology to
understand properties of the Universe, for it can thus give information on
both Gravity and neutrinos.

Studying neutrinos as quantum spinors in curved spacetime have been done in
several theoretical contexts, such as near Schwarzschild Blackholes [20], in
extended theories of gravity [21] or to derive fundamental uncertainty
relations [22, 23]. More specifically to dynamical DE, neutrinos have been
investigated in the context of mass-varying neutrinos [24, 25, 26, 27, 28],
pseudo-Dirac particles [29, 30] or Lorentz/$CPT$ violating theories [31, 32],
where $CPT$ stands for Charge Conjugation($C$), Parity($P$) and Time
reversal($T$) symmetries. It would be interesting therefore to expand on these
works to produce an observational trace of this kind of interactions.

As a first step along this way, one of us proposed a model, called DEν, in
which a scalar field is “frozen” in place via an interaction with neutrino
[33]. This model, by construction, mimics a cosmological constant from the
point of view of cosmological observables (expansion history and
perturbations) and thus it does not leave any significant imprint in these
classical observables. We then expanded upon that work [34], and looked at DEν
in the context of quantum spinors in curved spacetime, in addition to beyond
Standard Model scalar-spinor interactions.

In this work, we further develop ref. [34], as well as the works previously
mentioned, to produce an observable that could be measured experimentally,
which can then differentiate between DE models. We look at a more general
massive spinor-scalar field interaction(section 2), and then derive a
generalized formula for the oscillation probability in an arbitrary
spacetime(section 3). Although a scalar field DE scenario does not include all
types of DE, nevertheless it incorporates a large class of DE models,
including scalar-tensor theories of gravity such as Horndeski [35, 36]. That
is the reason why we focus on such an interaction here. Afterwards, we specify
to cosmological constant DE(section 4.1), what is known as $\Lambda$CDM model,
and quintessence [37, 38] with a neutrino-scalar interaction as presented in
the DEν model(section 4.2). For the latter, we look at various neutrino and
DEν properties and compare them to the former model. We finish by presenting a
summary and future prospects(section 5).

Throughout the paper, we use units in which $\hbar=c=1$, where $c$ is the
speed of light and $\hbar$ is the reduced Planck constant. Moreover, the
metric signature we use is the mostly positive one, $(-,+,+,+)$, and Greek
indices will be used for spacetime coordinates $(0,1,2,3)$, while Latin ones
are dedicated for spatial coordinates only, $(1,2,3)$. In addition, for
neutrino states notations, we use Greek and Latin indices to describe flavor
and mass states, respectively.

##  2 Generalized Formalism

Let us consider a general interaction of a spinor field, $\psi$, with a
classical scalar field, $\varphi$, in curved spacetime with metric
$g_{\mu\nu}$. This is described by the following action:

$\displaystyle S=\int d^{4}x\sqrt{-g}\bigg{[}$
$\displaystyle\frac{1}{2}R-\frac{1}{2}\mathcal{D}_{\mu}\mathcal{\varphi}{D}^{\mu}\varphi-V(\varphi)$
$\displaystyle+i\big{(}\bar{\psi}\gamma^{\mu}\mathcal{D}_{\mu}\psi-\mathcal{D}_{\mu}\bar{\psi}\gamma^{\mu}\psi\big{)}-2m\bar{\psi}\psi+\zeta\Theta\bigg{]}.$
(1)

In here, $g$ is the determinant of $g_{\mu\nu}$ and $R$ is the Ricci scalar,
the trace of the Ricci tensor $R_{\mu\nu}$. Moreover, $\mathcal{D}_{\mu}$ is a
generalized covariant derivative for fields with different spins in curved
spacetime (see Ref.[39, 40, 41] for more details on quantum fields in curved
backgrounds). For example, when acting a particle of spin 0,
$\mathcal{D}_{\mu}$ reduces to $\partial_{\mu}$, the usual partial derivative
in flat spacetime. We will see shortly the form it takes when acting on
spinors. Another term that appears in (1) is the scalar field potential
$V(\varphi)$, which describes self interactions of $\varphi$. Furthermore,
$\gamma^{\mu}=\\{\gamma^{0},\gamma^{1},\gamma^{2},\gamma^{3}\\}$ are the four
Dirac matrices, $\bar{\psi}=\psi^{\dagger}\gamma^{0}$, with $\psi^{\dagger}$
being the complex transpose of $\psi$, and $m$ is the spinor field’s mass.
Finally, $\zeta$ is the coupling constant for the general interaction term
between $\psi$ and $\varphi$,
$\Theta\big{(}\psi,\bar{\psi},\varphi,X^{\mu}_{\psi},X^{\mu}_{\bar{\psi}},X^{\mu}_{\varphi}\big{)}$,
with
$X^{\mu}_{\psi}\big{(}X_{\bar{\psi}}^{\mu}\big{)}=\mathcal{D}^{\mu}\psi\big{(}\mathcal{D}^{\mu}\bar{\psi}\big{)}$
and $X_{\varphi}^{\mu}=\partial^{\mu}\varphi$.

By setting the variational derivative of the action (1) with respect to(w.r.t)
$\bar{\psi}$ to 0, we get the modified Dirac equation for $\psi$, i.e.

$\frac{1}{\sqrt{-g}}\frac{\delta
S}{\delta\bar{\psi}}=0\quad\Rightarrow\quad\big{(}i\gamma^{\mu}\mathcal{D}_{\mu}-m\big{)}\psi=-\frac{\zeta}{2}\bigg{(}\frac{\partial\Theta}{\partial\bar{\psi}}-\mathcal{D}^{\mu}\frac{\partial\Theta}{\partial
X_{\bar{\psi}}^{\mu}}\bigg{)}\equiv-\frac{\zeta}{2}\frac{\delta\Theta}{\delta\bar{\psi}}.$
(2)

At this stage, one could say that the interaction is of the most general form,
for neither $\Theta$ nor the metric have been specified. However, as we are
considering a more phenomenological approach to the question at hand, it would
be more useful to look for a practical form of the variational derivative of
$\Theta$ w.r.t $\bar{\psi}$. This would allow us to calculate observables that
could be eventually measured by experiments. Moreover, it would prove useful
to divide the interaction term into flavor-invariant and flavor-dependent
parts, to study how each would affect the transition probability from one
flavor state to another. Intuitively, one would expect that the former should
not modify this flavor oscillation probability, since by definition the latter
is a transition between flavors. However, as we will see in the next section,
this is not always the case.

There are several ways in which one can implement these considerations. For
instance, in order to have the right hand side(r.h.s) of (2) mathematically
and dimensionally consistent with its left hand side(l.h.s), one possibility
is:

$-\frac{\zeta}{2}\frac{\delta\Theta}{\delta\bar{\psi}}=\big{(}\xi\gamma^{\mu}F_{\mu}(\varphi,X_{\varphi}^{\mu})+\xi_{f}G(\varphi,X_{\varphi}^{\mu})\big{)}\psi.$
(3)

The first term on the r.h.s is a global interaction in flavor space, i.e. it
couples to all flavors with the same strength $\xi$. On the other hand, the
second term is a flavor-specific term, with the coupling strength $\xi_{f}$
depending on which flavor is being considered. Another term that could be
added is a kinetic coupling, such as $\xi
F(\varphi,X_{\varphi}^{\mu})\gamma^{\mu}\mathcal{D}_{\mu}\psi$. However, as
such a term could produce effects on other Cosmological observables(see
appendix in [34]), we will not be considering it here. Finally, for the
purpose we seek of studying $\Lambda$CDM and quintessence(in the context of
DEν), it turns out that having

$-\frac{\zeta}{2}\frac{\delta\Theta}{\delta\bar{\psi}}=\big{(}\xi
F(\varphi,X_{\varphi}^{\mu})+\xi_{f}\gamma^{\mu}G_{\mu}(\varphi,X_{\varphi}^{\mu})\big{)}\psi.$
(4)

is more useful, and will therefore be used in the following sections

##  3 Neutrino Oscillation

In this section, for simplicity, we will be studying two-flavor neutrino
oscillation in curved space-time, although one could generalize the analysis
to three-flavor oscillations (see Ref. [20, 21, 22] for more details on
neutrino oscillations in curved space-time). In addition, a more stringent
study of neutrino oscillations in curved backgrounds would rely on the full
quantum field theoretic treatment (see Ref. [23, 42] and references within for
further details). However, for the purpose of studying neutrino interaction
with DE, the quantum mechanical treatment presented here is sufficient for
comparison with observations. We will look at the quantum field theoretic
treatment for both fields in future works.

###  3.1 Transition Amplitude’s Evolution

The first step in studying neutrino oscillations is to expand a state of
flavor $\alpha$, $|\nu_{\alpha}\rangle$, in terms of mass eigenstates,
$|\nu_{j}\rangle$:

$|\nu_{\alpha}(\lambda)\rangle=\sum_{j=1,2}U_{\alpha
j}e^{i\Phi(\lambda)}|\nu_{j}\rangle,$ (5)

where $\lambda$ is the monotonically increasing affine parameter along the
neutrino world-line and $U_{\alpha j}$ is the two-flavor mixing matrix, given
by:

$\displaystyle U$ $\displaystyle=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}$ (6)

with $\theta$ being the mixing angle. Moreover, $\Phi(\lambda)$ is the mass
eigenstate’s evolution operator [20]:

$\Phi(\lambda)=\int_{\lambda_{0}}^{\lambda}P_{\mu}\frac{dx^{\mu}}{d\lambda^{\prime}}d\lambda^{\prime}$
(7)

where $P^{\mu}$ is the 4-momentum operator, $dx^{\mu}/d\lambda$ is a null
vector tangent to the neutrino world-line, and $\lambda_{0}$($\lambda$) is the
affine parameter’s value at the observer(source).

One can see that eq. (5) is a solution for the Schrödinger like equation

$i\frac{d}{d\lambda}|\nu_{\alpha}(\lambda)\rangle=\Phi(\lambda)|\nu_{\alpha}(\lambda)\rangle,$
(8)

and therefore the transition amplitude between states $\alpha$ and $\beta$,

$\Psi_{\alpha\beta}=\langle\nu_{\beta}|\nu_{\alpha}(\lambda)\rangle,$ (9)

satisfies

$i\frac{d}{d\lambda}\Psi_{\alpha\beta}=\Phi(\lambda)\Psi_{\alpha\beta}.$ (10)

The ultimate goal is to find the transition probability between flavors
$\alpha$ and $\beta$, i.e.

$P_{\beta\rightarrow\alpha}=|\Psi_{\alpha\beta}|^{2}=|\langle\nu_{\beta}|\nu_{\alpha}(\lambda)\rangle|^{2},$
(11)

and therefore we need to calculate $\Phi(\lambda)$, or more specifically,
$P_{\mu}dx^{\mu}/d\lambda$, as has been pointed out before [20, 34].

For the purpose at hand, let us start with a system of two-flavors, electron
and muon neutrinos $\nu_{e}$ and $\nu_{\mu}$, respectively, that is
$\alpha,\beta=e,\mu$. Let

$\displaystyle\boldsymbol{\psi}$ $\displaystyle=\begin{pmatrix}\psi_{e}\\\
\psi_{\mu}\end{pmatrix},$ (12)

be a vector of spinor fields. The modified Dirac equation (2) for this system
becomes:

$\bigg{(}i\gamma^{\mu}\mathcal{D}_{\mu}-\mathcal{M}_{f}\bigg{)}\boldsymbol{\psi}=\big{(}\xi
F(\varphi,X_{\varphi}^{\mu})+\xi_{f}\gamma^{\mu}G_{\mu}(\varphi,X_{\varphi}^{\mu})\big{)}\boldsymbol{\psi}$
(13)

where the vacuum mass matrix in flavor space is given by

$\displaystyle\mathcal{M}_{f}^{2}$
$\displaystyle=U\begin{pmatrix}m_{1}^{2}&0\\\
0&m_{2}^{2}\end{pmatrix}U^{\dagger}$ (14)

with $m_{1},m_{2}$ being the masses of mass states
$|\nu_{1}\rangle,|\nu_{2}\rangle$, respectively.

Now we introduce the explicit form of the covariant derivative
$\mathcal{D}_{\mu}$. As explained in Ref. [39, 40], when studying spinors in
curved spacetime, one needs to introduce a local inertial coordinate system,
with its own set of Dirac matrices $\gamma^{a}$, and link it to the general
one using tetrad fields $e^{\mu}_{a}$, where Greek indices correspond to
general coordinates, while Latin ones for the local system. With this, we can
write

$\gamma^{\mu}\mathcal{D}_{\mu}=\gamma^{a}e^{\mu}_{a}\big{(}\partial_{\mu}+\Gamma_{\mu}\big{)}$
(15)

where

$\Gamma_{\mu}=\frac{1}{8}\big{[}\gamma^{b},\gamma^{c}\big{]}e_{b}^{\nu}\nabla_{\mu}e_{c\nu}$
(16)

is called the spin-connection which takes into account the gravitational
effect on the particle’s spin. In eq. (16), $[\gamma^{b},\gamma^{c}]$ is the
commutator of $\gamma^{a}$ and $\gamma^{b}$, and $\nabla_{\mu}$ is the usual
covariant derivative of General Relativity [43]. With these relations, one can
then show that

$\gamma^{a}e_{a}^{\mu}\Gamma_{\mu}=i\gamma^{a}e_{a}^{\mu}A_{G\mu}$ (17)

with

$A_{G}^{\mu}=\frac{1}{4}\sqrt{-g}e_{a}^{\mu}\epsilon^{abcd}(\partial_{\sigma}e_{b\nu}-\partial_{\nu}e_{b\sigma})e_{c}^{\nu}e_{d}^{\sigma}$
(18)

where $\epsilon^{abcd}$ is the local four dimensional Levi-Civita symbol.
Inserting eqs. (15), (17) and (18) in eq. (13) and moving all terms to the
l.h.s, we get:

$\bigg{\\{}i\gamma^{\mu}\bigg{[}\partial_{\mu}+i\bigg{(}A_{G\mu}+\xi_{f}G_{\mu}\bigg{)}\bigg{]}-\bigg{[}\mathcal{M}_{f}-\xi
F\bigg{]}\bigg{\\}}\boldsymbol{\psi}=0.$ (19)

In order to get non-trivial solutions for the above system, the determinant of
the braces must be 0. This results in a modified mass-shell relation:

$\bigg{(}P^{\mu}+A^{\mu}\bigg{)}\bigg{(}P_{\mu}+A_{\mu}\bigg{)}=\tilde{\mathcal{M}}_{f}^{2},$
(20)

where $A^{\mu}=A_{G}^{\mu}+\xi_{f}G_{\mu}$ and

$\displaystyle\tilde{\mathcal{M}}_{f}^{2}$
$\displaystyle=U\begin{pmatrix}\tilde{m}_{1}^{2}&0\\\
0&\tilde{m}_{2}^{2}\end{pmatrix}U^{\dagger}$ (21)

with $\tilde{m}_{i}=m_{i}-\xi F$ for $i=1,2$. It should be noted here that
$\tilde{m}$ is not a mass-varying neutrino, rather an effective mass due to
the interaction with another field (see [24, 25, 26] for comparison). From eq.
(20), one can show that

$P_{\mu}\frac{dx^{\mu}}{d\lambda}=\frac{1}{2}\tilde{\mathcal{M}}_{f}^{2}-\frac{dx^{\mu}}{d\lambda}A_{\mu},$
(22)

which finally implies, from eqs. (7) and (10), that

$i\frac{d}{d\lambda}\Psi_{\alpha\beta}=\bigg{[}\frac{1}{2}\tilde{\mathcal{M}}_{f}^{2}+V_{I}\bigg{]}\Psi_{\alpha\beta},$
(23)

with $V_{I}=-A_{\mu}dx^{\mu}/d\lambda$. In deriving eq. (22), two well
motivated assumptions have been made based on the fact that we are focusing on
high energy neutrinos [20, 34]. First, we consider neutrinos as energy
eigenstates, i.e. $P^{0}=dx^{0}/d\lambda$, and second, $P^{i}$ and
$dx^{i}/d\lambda$ are assumed parallel, that is
$P^{i}=(1-\varepsilon)dx^{i}/d\lambda$, with $\varepsilon\ll 1$ for high-
energy neutrinos111The second condition can be relaxed since we are eventually
taking the inner product of the two vectors, so that the perpendicular part
does not contribute..

###  3.2 Transition Probability

Let us now be more explicit, and look at each component of eq. (23). With some
matrix algebra, it can be shown that

$\displaystyle\tilde{\mathcal{M}}_{f}^{2}$
$\displaystyle=\bigg{(}\tilde{m}_{1}^{2}+\frac{1}{2}\tilde{\Delta}\bigg{)}I+\frac{1}{2}\tilde{\Delta}\begin{pmatrix}-\cos
2\theta&\sin 2\theta\\\ \sin 2\theta&\cos 2\theta\end{pmatrix},$ (24)

where $I$ is the $2\times 2$ identity matrix in flavor space and
$\tilde{\Delta}_{m}^{2}=\tilde{m}_{2}^{2}-\tilde{m}_{1}^{2}=\Delta_{m}^{2}-2\xi
F\Delta_{m}$(up to 1${}^{\text{st}}$ order in $\xi$), with
$\Delta_{m}^{2}=m_{2}^{2}-m_{1}^{2}$ and $\Delta_{m}=(m_{2}-m_{1})$. Note the
resemblance to the MSW effect [44, 45], with the difference being that
interactions with matter are substituted by those with spacetime and DE. Also,
it is safe to ignore the term proportional to $I$ in eq. (24), since it is
common for both transition amplitudes, and therefore will cancel when we
calculate the probability. If we start initially from a $\nu_{e}$ state, for
instance, the evolution equation for the transition amplitudes becomes:

$\displaystyle i\frac{d}{d\lambda}\begin{pmatrix}\Psi_{ee}\\\
\Psi_{e\mu}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}-\frac{1}{4}\tilde{\Delta}_{m}^{2}\cos
2\theta+\xi_{e}V_{I}&\frac{1}{4}\tilde{\Delta}_{m}^{2}\sin 2\theta\\\
\frac{1}{4}\tilde{\Delta}_{m}^{2}\sin
2\theta&\frac{1}{4}\tilde{\Delta}_{m}^{2}\cos
2\theta+\xi_{\mu}V_{I}\end{pmatrix}\begin{pmatrix}\Psi_{ee}\\\
\Psi_{e\mu}\end{pmatrix}$
$\displaystyle\equiv\boldsymbol{M}\begin{pmatrix}\Psi_{ee}\\\
\Psi_{e\mu}\end{pmatrix}.$ (25)

Notice that the gravitational contribution $A_{G\mu}$ has been dropped from
the interaction term. This is because it is proportional to $I$ in flavor
space, and therefore does not contribute to the oscillation probability [20].
In addition to that, in spatially homogeneous and isotropic universes, such as
FRW, this term is 0 identically [34].

From eq. (25), we can proceed by diagonalizing $\boldsymbol{M}$, which has

$v_{\pm}=\frac{1}{4}\bigg{[}2\big{(}\xi_{e}+\xi_{\mu}\big{)}V_{I}\pm\sqrt{\big{[}\tilde{\Delta}_{m}^{2}\cos
2\theta-2V_{I}(\xi_{e}-\xi_{\mu})\big{]}^{2}+(\tilde{\Delta}_{m}^{2})^{2}\sin^{2}2\theta}\bigg{]}$
(26)

as eigenvalues, and

$\displaystyle\tilde{U}$
$\displaystyle=\begin{pmatrix}\cos\tilde{\theta}&\sin\tilde{\theta}\\\
-\sin\tilde{\theta}&\cos\tilde{\theta}\end{pmatrix}$ (27)

as the unitary matrix that diagonalizes it, with222One way of deriving eq.
(28) is to perform the matrix product $\tilde{U}^{T}\boldsymbol{M}\tilde{U}$,
and equate it to $\text{diag}\\{v_{-},v_{+}\\}$.

$\displaystyle\cos 2\tilde{\theta}=\frac{\tilde{\Delta}_{m}^{2}\cos
2\theta-2V_{I}\big{(}\xi_{e}-\xi_{\mu}\big{)}}{\sqrt{\big{[}\tilde{\Delta}_{m}^{2}\cos
2\theta-2V_{I}\big{(}\xi_{e}-\xi_{\mu}\big{)}\big{]}^{2}+(\tilde{\Delta}_{m}^{2})^{2}\sin^{2}2\theta}};$
$\displaystyle\sin 2\tilde{\theta}=\frac{\tilde{\Delta}_{m}^{2}\sin
2\theta}{\sqrt{\big{[}\tilde{\Delta}_{m}^{2}\cos
2\theta-2V_{I}\big{(}\xi_{e}-\xi_{\mu}\big{)}\big{]}^{2}+(\tilde{\Delta}_{m}^{2})^{2}\sin^{2}2\theta}}.$
(28)

In analogy with the flavor-mass bases transformation, let us define

$\displaystyle\boldsymbol{\phi}_{e}$
$\displaystyle\equiv\begin{pmatrix}\phi_{e-}\\\
\phi_{e+}\end{pmatrix}=\tilde{U}^{T}\begin{pmatrix}\Psi_{ee}\\\
\Psi_{e\mu}\end{pmatrix},$ (29)

as a vector of transition amplitudes in an effective mass basis,
$\\{\nu_{-},\nu_{+}\\}$, that takes into account the neutrino interaction with
gravity and DE. Using the unitarity of $\tilde{U}$ and eq. (25), it can be
shown that $\boldsymbol{\phi}_{e}$ satisfies:

$\displaystyle i\frac{d}{d\lambda}\begin{pmatrix}\phi_{e-}\\\
\phi_{e+}\end{pmatrix}=\begin{pmatrix}v_{-}&-i\frac{d\tilde{\theta}}{d\lambda}\\\
i\frac{d\tilde{\theta}}{d\lambda}&v_{+}\end{pmatrix}\begin{pmatrix}\phi_{e-}\\\
\phi_{e+}\end{pmatrix}.$ (30)

Notice that the off-diagonal terms come from transforming the l.h.s of eq.
(25). Also, in the case where there is no mixing between effective mass
states, then $d\tilde{\theta}/d\lambda=0$, and the transition amplitudes
evolve as:

$\phi_{ej}=\big{(}\cos\omega_{j}+i\sin\omega_{j}\big{)}\phi_{ej}(0)$ (31)

for $j=+,-$, where $\phi_{ej}(0)$ is the initial condition and

$\omega_{j}(\lambda)=\int_{\lambda_{0}}^{\lambda}v_{j}d\lambda^{\prime}.$ (32)

This is known as the adiabatic evolution condition which, as we will see in
the next section, applies to the DE scenarios we will examine. Further, one
important consequence of adiabaticity is that the flavor-specific interaction
will be constant along the neutrino’s world-line. To see this, differentiate
$\cos 2\tilde{\theta}$ from eq. (28) w.r.t $\lambda$:

$\frac{d\tilde{\theta}}{d\lambda}=\frac{\sin
2\tilde{\theta}}{\Delta_{m}}\biggl{\\{}-2\xi\Delta_{m}\frac{dF}{d\lambda}\bigg{[}\frac{\cos
2\tilde{\theta}\sin 2\theta}{\sin 2\tilde{\theta}}-\cos
2\theta\bigg{]}+\frac{dV_{I}}{d\lambda}(\xi_{e}-\xi_{\mu})\biggr{\\}}.$ (33)

Since we expect gravitational and DE effects to be small compared to the
vacuum oscillations, that is $\xi,\xi_{f}\ll 1$, we can keep terms up to first
order in these coupling constants. With this assumption, by setting eq. (33)
to 0, we find that $dV_{I}/d\lambda=0$. This is again another analogy with the
MSW effect, where in adiabatic oscillations the interaction term is constant
along the path [46].

The final ingredient we need to get the oscillation probability is initial
conditions. As we are considering an initial $\nu_{e}$ state, we can write

$\displaystyle\boldsymbol{\phi}_{e}(0)$
$\displaystyle\equiv\begin{pmatrix}\phi_{e-}(0)\\\
\phi_{e+}(0)\end{pmatrix}=\begin{pmatrix}\cos\tilde{\theta}&-\sin\tilde{\theta}\\\
\sin\tilde{\theta}&\cos\tilde{\theta}\end{pmatrix}\begin{pmatrix}1\\\
0\end{pmatrix}=\begin{pmatrix}\cos\tilde{\theta}\\\
\sin\tilde{\theta}\end{pmatrix},$ (34)

where, by construction, an initial $\nu_{e}$ state corresponds to
$\Psi_{ee}(0)=1$ and $\Psi_{e\mu}(0)=0$. By acting with the inverse
transformation of eq. (29), we can calculate the amplitude $\Psi_{e\mu}$, and
thus, with the initial conditions eq. (34), we finally obtain the
$\nu_{e}\rightarrow\nu_{\mu}$ transition probability:

$P_{\nu_{e}\rightarrow\nu_{\mu}}=|\Psi_{e\mu}|^{2}=\sin^{2}2\tilde{\theta}\sin^{2}\bigg{(}\frac{\omega_{-}-\omega_{+}}{2}\bigg{)}.$
(35)

Let us now look in more detail into the oscillating term in eq. (35). If we
use the above mentioned approximation ($\xi$, $\xi_{f}\ll 1$), one can show
that

$\displaystyle\omega_{-}-\omega_{+}$
$\displaystyle=\int_{\lambda_{0}}^{\lambda}(v_{-}-v_{+})d\lambda^{\prime}$
$\displaystyle\approx\frac{\Delta_{m}^{2}}{2}(\lambda_{0}-\lambda)+V_{I}\cos
2\theta(\xi_{e}-\xi_{\mu})(\lambda-\lambda_{0})+\xi\Delta_{m}\int_{\lambda_{0}}^{\lambda}Fd\lambda^{\prime}.$
(36)

The first term in eq. (36) corresponds to the usual vacuum oscillation term.
Indeed, if one neglects the interactions completely, i.e. $F=G=0$, and
consider Minkowski spacetime, that is $d\lambda=dt/E=dx/E$, we get

$\omega_{-}-\omega_{+}=\omega_{\text{std}}\equiv\frac{\Delta_{m}^{2}L}{2E}$
(37)

where $L$ is the distance traveled by the neutrinos. This results in

$P^{\text{std}}_{\nu_{e}\rightarrow\nu_{\mu}}=\sin^{2}2\theta\sin^{2}\bigg{(}\frac{\Delta_{m}^{2}L}{4E}\bigg{)},$
(38)

which is the standard vacuum transition probability in flat spacetime [46].
The second term in eq. (36) is the flavor-specific correction, and the third
term is an integrated correction from the flavor-invariant interaction. This
is where we see that the latter does affect the transition probability, both
in amplitude, through $\sin 2\tilde{\theta}$ in eq. (35), and in period.

Let us finish the analysis by writing a Signal-to-Noise-like expression for
the oscillation probability eq. (35). If we substitute the expression for
$\sin 2\tilde{\theta}$ from eq. (28) into eq. (35), and then expand all
functions of the interactions up to 1${}^{\text{st}}$ order, we get

$\frac{\delta
P}{P}\equiv\frac{P_{\nu_{e}\rightarrow\nu_{\mu}}-P^{\text{vac}}_{\nu_{e}\rightarrow\nu_{\mu}}}{P^{\text{vac}}_{\nu_{e}\rightarrow\nu_{\mu}}}=\frac{4V_{I}(\xi_{e}-\xi_{\mu})}{\Delta_{m}^{2}}+\cot\bigg{(}\frac{1}{2}\omega_{\text{vac}}\bigg{)}\omega_{\text{DE}},$
(39)

where

$P^{\text{vac}}_{\nu_{e}\rightarrow\nu_{\mu}}=\sin^{2}2\theta\sin^{2}\omega_{\text{vac}}$
(40)

is the transition probability in vacuum, with frequency

$\omega_{\text{vac}}=\frac{\Delta_{m}^{2}}{2}(\lambda_{0}-\lambda),$ (41)

and

$\omega_{\text{DE}}=V_{I}\cos
2\theta(\xi_{e}-\xi_{\mu})(\lambda-\lambda_{0})+\xi\Delta_{m}\int_{\lambda_{0}}^{\lambda}Fd\lambda^{\prime}$
(42)

is the additional contribution to the oscillation frequency due to the
interaction with DE.

In this section, we looked at how a type of general interactions between
neutrinos and DE, in a generic spacetime, can affect the probability of
oscillations, with the final result given in eq. (35). Now we can specify the
interaction to known DE models, particularly a cosmological constant and
scalar field based DE, and thus establish the distinction between them.

##  4 Oscillation Probability for Specific DE Models

As mentioned in the Introduction 1, we have focused on the interaction of
neutrinos with scalar fields since the latter includes a large class of DE
models, such as some modified gravity scenarios and Quintessence. Having
established a general formalism for the interaction of neutrinos with a scalar
field in the previous section, we will now focus on two DE energy models: a
Cosmological Constant $\Lambda$ and Quintessence.

###  4.1 $\Lambda$CDM

This model is the simplest model describing our Universe, and has sustained a
great deal of observational test [3, 47]. Taking GR as the theory of gravity,
$\Lambda$CDM has two main components in the late universe: a cosmological
constant DE, $\Lambda$, and cold Dark Matter(CDM). The metric of spacetime
that best describes it is FLRW:

$ds^{2}=-dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j},$ (43)

where $t$ is cosmic time, $x^{i}$, for $i=1,2,3$, are comoving spatial
coordinates, $\delta_{ij}$ is the Kronecker delta and $a(t)$ is the scale
factor that incorporates the universe’s expansion. The resulting evolution
equation for the scale factor will be the usual first Friedmann equation:

$H^{2}=H_{0}^{2}\big{(}\Omega_{m_{0}}a^{-3}+\Omega_{\Lambda_{0}}\big{)},$ (44)

where $H=\dot{a}/a$ is the Hubble parameter, with $H_{0}$ being its value
today and $\Omega_{m0}$, $\Omega_{\Lambda 0}$ are today’s matter and $\Lambda$
density parameters, respectively, with the most recent measurement given by
the Planck Collaboration [3]. Note that this equation is not altered for two
reasons. First, we are considering gravity in terms of the background
spacetime, and not from a quantum perspective, hence Einstein equations, from
which eq. (44) is derived, are still the same. Second, the neutrino density
parameter has not been added since it is small compared to that of matter and
$\Lambda$ [48].

To see the effect of DE in this model on neutrino oscillations, let us start
by noting that when DE is $\Lambda$, $\xi=\xi_{f}=0$ in eq. (4), implying the
automatic satisfaction of the adiabaticity condition,
$d\tilde{\theta}/d\lambda=0$, from eq. (33), and thus
$\sin\tilde{\theta}(\cos\tilde{\theta})=\sin\theta(\cos\theta)$. Therefore,
from the first equality in eq. (36), we define

$\omega_{\Lambda}\equiv\omega_{-}-\omega_{+}\big{|}_{DE=\Lambda}=\frac{\Delta_{m}^{2}}{2}\int_{t_{\text{em}}}^{t_{0}}\frac{1}{E}dt,$
(45)

with the second equality meaning eq. (36) when DE is $\Lambda$. Here, $t_{0}$
is today, $t_{\text{em}}$ is the time of neutrino emission and $E=dt/d\lambda$
is the 0${}^{\text{th}}$ component of the null tangent vector
$dx^{\mu}/d\lambda$, which is also the neutrino’s energy. Since the latter
follows the geodesic equation, as shown in [20, 34], then $E=E_{0}/a$, with
$E_{0}$ being the neutrino energy at detection, and thus, using eq. (44),

$\displaystyle\omega_{\Lambda}$
$\displaystyle=\frac{\Delta_{m}^{2}}{2H_{0}E_{0}}\int_{a_{\text{em}}}^{1}\bigg{(}\Omega_{m_{0}}a^{-3}+\Omega_{\Lambda_{0}}\bigg{)}^{-1/2}da$
$\displaystyle=\frac{\Delta_{m}^{2}}{2H_{0}E_{0}}\int_{0}^{z_{\text{em}}}\bigg{(}\Omega_{m_{0}}(1+z)^{7}+\Omega_{\Lambda_{0}}(1+z)^{4}\bigg{)}^{-1/2}dz,$
(46)

where $a_{\text{em}}(z_{\text{em}})$ is the scale factor(redshift) at neutrino
emission, and with the usual normalization $a_{0}=1$and $z_{0}=0$. On the
other hand, if one takes a simple approach (SA) to neutrino oscillations in an
expanding universe, and substitutes $L$ and $E$ in eq. (38) by the luminosity
distance,

$D_{L}=(1+z_{e})H_{0}^{-1}\int_{0}^{z_{e}}\bigg{(}\Omega_{m_{0}}(1+z)^{3}+\Omega_{\Lambda_{0}}\bigg{)}^{-1/2}dz$
(47)

and $E=E_{0}(1+z_{e})$, respectively, we get,

$\omega_{\text{SA}}=\frac{\Delta_{m}^{2}}{2H_{0}E_{0}}\int_{0}^{z_{e}}\bigg{(}\Omega_{m_{0}}(1+z)^{3}+\Omega_{\Lambda_{0}}\bigg{)}^{-1/2}dz.$
(48)

Finally, inserting eqs. (46) and (48) in eq. (35) gives the two-flavor
neutrino oscillation probability in the $\Lambda$CDM model,

$P_{\Lambda}=\sin^{2}2\theta\sin^{2}\omega_{\Lambda},$ (49)

and in the SA,

$P_{\text{SA}}=\sin^{2}2\theta\sin^{2}\omega_{\text{SA}}.$ (50)

In Figure 1, we plot the evolution of $P_{\Lambda}$(solid black curve) and
$P_{\text{SA}}$(dotted blue curve) as a function of redshift, to compare the
two approaches. To this end, we took $\Delta_{m}^{2}=7.53\times
10^{-5}$eV2333Here we used mass states 1 and 2 from Ref. [49] as ours. One can
check that other values of $\Delta_{m}^{2}$ reported there does not alter the
evolution of the frequencies eqs.(46, 48). Physically, this is due to the
absence of a direct interaction between DE and neutrinos. Mathematically, this
is because the coefficient multiplying the integrals in eqs.(46, 48) includes
$H_{0}^{-1}\sim\mathcal{O}(10^{33})$eV, which wipes out the
$\mathcal{O}(10^{2})$eV2 difference between $\Delta_{m}^{2}$s. and
$E_{0}=10^{16}$eV, a value to which neutrino detectors are on average
sensitive to [50]. Further, we used
$\Omega_{m_{0}}=0.315,\Omega_{\Lambda_{0}}=0.685$ and $H_{0}=1.44\times
10^{-33}$eV as reported in [3]. For redshifts higher than $\sim 2$, the
difference between the two probabilities stabilizes at around 80%, as can be
seen from figure 2. On the other hand, the latter shows, for the
observationally more interesting range of redshifts ($0\leq z\leq 0.5$), the
difference can reach up to 50% while they coincide for redshift 0, as
expected.

The difference between eqs.(46, 49) and eqs.(48,50) is being highlighted here
to insure that, when doing neutrino observations, one cannot directly
substitute $D_{L}(z)$ and $E(z)$ as neutrino traveling length and energy,
respectively. This will not properly take into account the evolution of a spin
1/2 particle in a curved background. Rather, one should use the formalism
presented in section 3, for a more general interaction with a scalar field in
curved spacetime, or eqs.(46, 49) for $\Lambda$CDM. The same idea applies to
other models of DE, however there will be differences in the evolution of the
oscillation probability, as we will see next.

Figure 1: The two-flavor neutrino oscillation probability, divided by
$\sin^{2}2\theta$, as a function of redshift of emission, $z_{e}\in[0,10]$,
for the two cases $P_{\Lambda}$(solid black curve) and $P_{\text{SA}}$(dotted
blue line), given by eqs. (45),(49) and eqs. (48),(50), respectively. To be
specific, due to the large value of $C=\Delta_{m}^{2}/(2H_{0}E_{0})$, we used
$\omega\mod{2\pi C}$ as the argument of $\sin^{2}$ in eqs.(49)-(50) to avoid
numerical instabilities. The values of the different parameters used in these
equations is given in the text below eq.(50). Figure 2: Absolute value of the
fractional difference between the oscillation probabilities $P_{\Lambda}$ and
$P_{\text{SA}}$, given by eqs. (45),(49) and eqs. (48),(50), respectively, as
a function of emission redshift, $z_{e}\in[0,10]$. To be specific, due to the
large value of $C=\Delta_{m}^{2}/(2H_{0}E_{0})$, we used $\omega\mod{2\pi C}$
as the argument of $\sin^{2}$ in eqs.(49)-(50) to avoid numerical
instabilities. Moreover, we shifted both $P_{\Lambda}$ and $P_{\text{SA}}$ by
$10^{-5}$ to avoid the singularity when they are 0 at $z_{e}=0$. The values of
different parameters used in these equations is given in the text below
eq.(50).

###  4.2 Quintessence

As a homogeneous canonical scalar field minimally coupled to gravity, Quint-
essence could be an explanation to the late time accelerated expansion [37,
51, 52, 53, 8]. One of the main reasons for introducing quintessence as an
alternative to a cosmological constant is to make DE dynamical, thereby
avoiding the cosmological constant and coincidence problems(see [38] and
references therein for more information on Quintessence).

In order to probe this model using neutrino oscillations, a coupling between
the scalar and spinor fields has to be introduced, otherwise the difference in
effect of quintessence and $\Lambda$ on the oscillation probability will be
difficult to observe. We consider the coupling introduced in [33], given the
name DEν model, and which we analyzed in [34]. In the present formalism, DEν
translates to $F=0$ and $G_{\mu}=\partial_{\mu}\varphi$ in eq. (4). As
mentioned in [33], such a derivative coupling is a low energy limit of the
model presented in [54], with the scalar field being a Nambu-Goldstone boson
resulting from the spontaneous symmetry breaking of Lepton number symmetry
[55, 56, 57]. This shows that such a coupling is motivated both from Particle
Physics and Cosmology points of view, hence it is being further analyzed here.

To start the analysis, recall that since the scalar field is homogeneous, its
energy density would be

$\rho_{\varphi}=\frac{1}{2}\dot{\varphi}^{2}+V(\varphi)$ (51)

where $\dot{\varphi}=\partial_{t}\varphi$ and $V_{\varphi}(\varphi)$ is the
potential energy of $\varphi$. Therefore, eq. (44) becomes

$H^{2}=H_{0}^{2}\big{(}\Omega_{m_{0}}a^{-3}+\Omega_{\varphi}\big{)},$ (52)

where

$\Omega_{\varphi}=\frac{8\pi G}{3H_{0}^{2}}\rho_{\varphi}$ (53)

is the density parameter of quintessence. Moreover, as already shown in [34],
this type of interactions does not affect the Klein-Gordon equation, which can
be written as:

$\frac{d}{da}\big{(}a^{6}\dot{\varphi}^{2}\big{)}=2a^{6}\frac{dV}{da}.$ (54)

For $\varphi$ to produce an accelerated expansion, it should satisfy the
condition:

$\dot{\varphi}^{2}\ll V(\varphi)\approx\rho_{\Lambda_{0}},$ (55)

where $\rho_{\Lambda_{0}}$ is the energy density of a cosmological constant
today. This means that, first, in eq. (52),
$\Omega_{\varphi}\approx\Omega_{\Lambda_{0}}$, and second, we can write444The
$7/2$ factor is to reduce numerical factors clustering.

$\frac{dV}{da}\equiv\frac{7}{2}\epsilon,$ (56)

where $\epsilon<\rho_{\Lambda_{0}}\sim\mathcal{O}(10^{-11}\text{eV}^{4})$555
Note that $\epsilon$ is not exactly the slow-roll parameter
$\varepsilon=d(H^{-1})/dt$, but one can show that $\varepsilon\approx
3\epsilon a/V$., and thus, from eq. (54), we get

$\dot{\varphi}=\sqrt{\epsilon a}.$ (57)

On the neutrino’s side, this type of interaction results in

$V_{I}=-\frac{dx^{\mu}}{d\lambda}G_{\mu}=-E\dot{\varphi}=-E_{0}\sqrt{\epsilon(1+z)},$
(58)

from which one can show, using eqs. (26), (28) and (33), that

$v_{\pm}=\frac{1}{2}V_{I}\big{[}\xi_{e}+\xi_{\mu}\mp\cos
2\theta(\xi_{e}-\xi_{\mu})\big{]}\pm\Delta_{m}^{2},$ (59) $\sin
2\tilde{\theta}=\sin
2\theta\bigg{[}1+\frac{4E_{0}\sqrt{\epsilon(1+z)}}{\Delta_{m}^{2}}\cos
2\theta(\xi_{e}-\xi_{\mu})\bigg{]}^{-1/2}$ (60)

and

$\frac{d\tilde{\theta}}{d\lambda}=\frac{\sin
2\tilde{\theta}}{\Delta_{m}}(\xi_{e}-\xi_{\mu})\frac{dV_{I}}{d\lambda},$ (61)

respectively. To check if the adiabaticity condition is satisfied for the
current case, differentiate eq. (58) w.r.t $\lambda$ and insert it in eq.
(61), to find

$\frac{d\tilde{\theta}}{d\lambda}\approx\frac{\sin
2\theta(\xi_{e}-\xi_{\mu})}{2\Delta_{m}}E_{0}^{2}\epsilon^{1/2}H_{0}\sqrt{\Omega_{m_{0}}(1+z)^{6}+\Omega_{\varphi_{0}}(1+z)^{3}}.$
(62)

From the fact that $H_{0}\sim\mathcal{O}(10^{-33})$eV [48],
$E_{0}\sim\mathcal{O}(10^{16})$eV(typical value for high-energy neutrinos
[50]), $\xi_{e,\mu}\sim\mathcal{O}(10^{-14})$eV-1 [34] and
$\epsilon\sim\mathcal{O}(10^{-11})$eV4, one can see that
$d\tilde{\theta}/d\lambda\ll v_{\pm}$, and therefore the adiabaticity
condition still holds, resulting in an oscillation frequency

$\displaystyle\omega_{Q}\approx$
$\displaystyle\frac{\Delta_{m}^{2}}{2E_{0}H_{0}}\int_{0}^{z_{e}}\sqrt{\frac{1+\frac{4E_{0}\sqrt{\epsilon(1+z)}}{\Delta_{m}^{2}}\cos
2\theta(\xi_{e}-\xi_{\mu})}{\Omega_{m_{0}}(1+z)^{7}+\Omega_{\varphi_{0}}(1+z)^{4}}}dz.$
(63)

Finally, from eq. (35), the two-flavor oscillation probability in the case
when DE is quintessence is:

$P_{Q}=\frac{\sin^{2}2\theta}{1+\frac{4E_{0}\sqrt{\epsilon(1+z)}}{\Delta_{m}^{2}}\cos
2\theta(\xi_{e}-\xi_{\mu})}\sin^{2}(\omega_{Q}/2).$ (64)

To study the difference between this model and $\Lambda$CDM, we plot(figure 3)
eq.(64) for values of $\xi_{i}$, $i=e,\mu$, ranging from $10^{-17}$ to
$10^{-14}$eV-1, in addition to eq.(49) for the $\Lambda$CDM case. We have
checked that smaller values of $\xi_{i}$ do not produce any noticeable
deviation from $\Lambda$CDM, while already at $10^{-14}$ we can see from
figure 3 that the deviation is $\sim 50\%$ at $z=2$. That is the reason why we
focus on this range of values of the couplings $\xi_{i}$. Moreover, we use the
same parameters used to produce figures 1 and 2(see text after eq. (50)), in
addition to $\cos 2\theta=0.4$ [49]. As the strength of the coupling
increases, the difference between the two models starts to become apparent at
redshift $\sim 0.5$, which is expected since then DE is becoming more
dynamical than in the case of $\Lambda$CDM.

In order to make the distinction between the different DE scenarios more
concrete, we study in more detail the dependence of $P_{Q}$ in eq. (64) on the
parameters $\epsilon$ and $\xi_{f}$ of this particular DE model. First, if
quintessence is slow-rolling, but not ultra slow-rolling, then
$\varepsilon$(see footnote 5) cannot be too small [58]. Taking
$\epsilon\in[10^{-14},10^{-12}]$, which corresponds to
$\varepsilon\in[10^{-3},10^{-1}]$, we find that
$\xi_{f}\in[10^{-15},10^{-13}]$ gives distinguishable stable results. On the
other hand, values beyond this interval would lead to unstable transition
probabilities. This shows that, in principle, the presented method here could
provide a complementary theoretical constraint to this type of coupling.

Second, if we now fix $\epsilon\sim\mathcal{O}(10^{-13})$, for instance, we
find that the difference between quintessence and $\Lambda$CDM starts to
become appreciable(i.e. more that a few %) for
$\xi_{f}\sim\mathcal{O}(10^{-14}-10^{-13})$. This is still consistent with
Particle Physics constraints for this type of coupling, which is
$\xi_{f}\lesssim 10^{-7}$eV-1 [59, 60]. Moreover, for values of $\xi_{f}$ that
differ from each other by at least half an order of magnitude, the transition
probabilities start deviating from each other by more than a few %. One may
conclude from this that there is a small window for fine-tuning in this model,
but not a too small one.

On another note, we can also explore how the results might change for
different neutrino parameters. Unlike for $\Lambda$CDM, the neutrino-
quintessence interaction is affected by the value of $\Delta_{m}^{2}$ and its
hierarchy, which is evident from the denominator of eq.(64). To see this, we
plot in the upper panel of figure 4 the probabilities shown in figure 3, but
for $\Delta_{m}^{2}=2.51\times 10^{-3}$eV2 (normal hierarchy), while in the
lower panel we use $\Delta_{m}^{2}=-2.56\times 10^{-3}$eV2 (inverted
hierarchy), with $\cos 2\theta\sim 0.2$ for both. These values of
$\Delta_{m}^{2}$ correspond to the difference between neutrino mass states 3
and 2 of the standard neutrino oscillation treatment [49].

There are a few things to be noted from these plots. First, if our neutrino
mass states are 3 and 2 from [49], it is more difficult to distinguish
$\Lambda$CDM from the quintessence model considered here for coupling
constants smaller than $10^{-14}$. This difficulty can be evaded once we
consider the full 3-flavor neutrino oscillations, which will be done in future
works. Our purpose here is merely to show that different DE models affect
neutrino oscillations differently. Second, even when we include all three
neutrino flavors, there will be a noticeable difference between the two
hierarchies for larger values of the coupling($\sim 10^{-14}$), as is apparent
from the two panels of figure 4. Therefore, such a neutrino-quintessence
interaction could require some fine-tuning to match future observations, which
puts it at equal, or less, footing with $\Lambda$CDM666Unless the value of the
coupling constant is derived from a more fundamental theory, which establishes
a fixed distinction between this model and $\Lambda$CDM..

Figure 3: Evolution of the neutrino oscillation probability with redshift in
the case of $\Lambda$CDM(solid black line) and quintessence, for neutrino-
quintessence coupling $\sim\mathcal{O}(10^{-17})$(dashed blue),
$10^{-16}$(dotted red), $10^{-15}$(dot-dash green) and $10^{-14}$(solid purple
line). The parameters used are given in the text, after eq. (50). Figure 4:
Evolution of the neutrino oscillation probability with redshift in the case of
$\Lambda$CDM(solid black line) and quintessence, for neutrino-quintessence
coupling $\sim\mathcal{O}(10^{-17})$(dashed blue), $10^{-16}$(dotted red),
$10^{-15}$(dot-dash green) and $10^{-14}$(solid purple line). The parameters
used are given in the text after eq. (50), except for $\cos 2\theta\sim 0.2$
and $\Delta_{m}^{2}$, which is $2.51\times 10^{-3}$eV2(upper panel) for normal
hierarchy, and $-2.56\times 10^{-3}$eV2(lower panel) for the inverted one.

###  4.3 Observational Strategy

Let us now comment on the relationship between our findings and observable
quantities. Note that, due to the fact that we are considering a two-flavor
neutrino system, direct comparison with neutrino observations would not be
very beneficial. Nevertheless, our main results, presented in figures 1 to 4,
do affect neutrino observations, and we will be exploring this in more detail
for three-flavor neutrinos in future work.

The main quantities observed at neutrino observatories, such as IceCube [50],
are neutrino fluxes. For instance, the electron neutrino flux, $F_{\nu_{e}}$
can be expressed as [61, 62]:

$F_{\nu_{e}}=\sum_{\alpha=e,\mu\tau}P_{\nu_{\alpha}\rightarrow\nu_{e}}F^{0}_{\nu_{\alpha}},$
(65)

where $F^{0}_{\nu_{\alpha}}$ is the flux of neutrinos with flavor $\alpha$ at
the source. It is in this expression that our results could affect neutrino
observations. The interaction of spinor neutrinos with curved spacetime will
alter this expression through the transition probability
$P_{\nu_{\alpha}\rightarrow\nu_{e}}$. More specifically to our case, depending
on which DE model is considered, eqs. (49) and (64) will give different
$P_{\nu_{\alpha}\rightarrow\nu_{e}}$ as a function of redshift, and thus the
neutrino flux detected will be different. Therefore, by calculating the
neutrino flux for each DE model, and compare it with observations, one can
distinguish between these models.

Another observational aspect worth mentioning is the experimental sensitivity
available for such effects to be observed. We would like first to highlight
that, when analyzing neutrino data, the usual emphasis is on the probability
and flux’s dependence on the neutrino’s energy. However, in addition to this
dependence, we are drawing attention here to the non-trivial effect of
spacetime curvature on the observational results, which in an FLRW context
translates into the dependence on the redshift. That is why in the analysis
above a value for the energy of$\sim$10PeV has been chosen. Such a value is
within reach of next generation neutrino detector IceCube-Gen2, which will
have a 5 times better sensitivity than IceCube [63].

##  5 Conclusions

In this paper, we have provided a proof of concept that distinct DE models can
be distinguished using neutrino oscillations, particularly through the
evolution of the oscillation probability with redshift. We first looked at a
more general interaction between two-flavor neutrinos, as quantum spinor
fields, and a classical scalar field in general spacetime. We focused on the
interaction with a scalar field since it comprises a large class of models for
DE, including $\Lambda$CDM, quintessence and some modified gravity scenarios
(such as Horndeski theory [35, 36]). Moreover, the interaction term considered
includes a part that couples equally to both neutrino flavors, a flavor-global
interaction, and another which is flavor dependent (see eq. (4)). The purpose
is to examine the different effect these two terms have on the oscillation
probability, which can be seen from the main result eqs. (35, 36) in section
3.2.

Furthermore, we applied this general formalism to two specific DE models,
$\Lambda$CDM and quintessence, to produce observable contrast between them
using neutrino oscillations. In the former model, we showed in figure 1 the
evolution of the oscillation probability with redshift when DE is a
cosmological constant. We also show in that figure the oscillation probability
in case of a direct substitution of the cosmological distance traveled by
neutrinos(such as the luminosity distance) and their energy in the standard
formula for neutrino oscillations eq (38), what we called SA. The point of
this contrast is that, if we detect $\nu_{e}$ from a type Ia supernova(SN),
for example, and we want to calculate their flux (which depends on the
$\nu_{e}$’s survival probability), SA would give a result $\sim
50-80\%$(depending on the SN’s redshift, see figure 2) more than the actual
value. This should be taken into account when doing neutrino observations in
the future [63, 64].

On the other hand, for quintessence, we looked at a derivative coupling
between neutrinos and the scalar field that is motivated by symmetry breaking
arguments [55, 56, 57], which was referred to in [33] as the DEν model. This
coupling, and others, have been already studied in [34], but we focused in
this work on the observational consequences of such a coupling which, without
it, $\Lambda$CDM and quintessence would be indistinguishable. In figure 3, we
show the oscillation probability’s evolution with redshift for the two models,
with the DEν coupling varied from $10^{-17}$ to $10^{-14}$eV-1. We also
investigated the effect several $\Delta_{m}^{2}$ values from Particle Physics
have on the probabilities, which in figure 3 was produced assuming mass states
1 and 2 from [49] as ours. This plot shows a clear distinction between
$\Lambda$CDM and quintessence for several values of the DEν coupling. However,
if we consider states 2 and 3 as our mass states, it would become more
difficult to distinguish the two DE models, unless the DEν coupling is at
least $\mathcal{O}(10^{-14})$, as seen in figure 4. Nevertheless, one can see
from the two plots of figure 4 a difference in the probability’s evolution
between the normal and inverted hierarchies, specially for high values of the
DEν coupling.

In the future, we would like to generalize the present work further, by
looking at the full three-flavor neutrino scenario, which should alleviate the
distinction between mass states choice previously mentioned. However, we
expect the difference between hierarchies’ choice to remain even in this case,
which prompts investigating its possible degeneracy with parameters of DE
models. Furthermore, one could also look at another type of general
interaction that could include other modified gravity models for DE, such as
extended gravity [65] or higher dimensions [66]. Finally, with the advancement
in neutrino detection techniques, we would expect these signals to appear in
near future terrestrial experiments [63, 64], or perhaps underground lunar
ones.

##  6 Acknowledgment

We would like to thank David Valcin and Nicola Bellomo for very helpful
discussions. The work of ARK and RJ is supported by MINECO grant
PGC2018-098866-B-I00 FEDER, UE. ARK and RJ acknowledge “Center of Excellence
Maria de Maeztu 2020-2023" award to the ICCUB (CEX2019- 000918-M).

## References

  * [1] S. Perlmutter, et al., Measurements of Omega and Lambda from 42 high redshift supernovae, Astrophys. J. 517 (1999) 565–586.  arXiv:astro-ph/9812133, doi:10.1086/307221.
  * [2] A. G. Riess, et al., Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116 (1998) 1009–1038.  arXiv:astro-ph/9805201, doi:10.1086/300499.
  * [3] N. Aghanim, et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6.  arXiv:1807.06209, doi:10.1051/0004-6361/201833910.
  * [4] The Dark Energy Spectroscopic Instrument (DESI), in: Bulletin of the American Astronomical Society, Vol. 51, 2019, p. 57.  arXiv:1907.10688.
  * [5] R. Laureijs, G. Racca, L. Stagnaro, J.-C. Salvignol, J. L. Alvarez, G. S. Criado, L. G. Venancio, A. Short, P. Strada, C. Colombo, G. Buenadicha, J. Hoar, R. Kohley, R. Vavrek, Y. Mellier, M. Berthe, J. Amiaux, M. Cropper, S. Niemi, S. Pottinger, A. Ealet, K. Jahnke, T. Maciaszek, F. Pasian, M. Sauvage, S. Wachter, U. Israelsson, W. Holmes, M. Seiffert, V. Cazaubiel, A. Anselmi, P. Musi, Euclid mission status, in: J. M. O. Jr., M. Clampin, G. G. Fazio, H. A. MacEwen (Eds.), Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, Vol. 9143, International Society for Optics and Photonics, SPIE, 2014, pp. 132 – 139.  doi:10.1117/12.2054883.
URL https://doi.org/10.1117/12.2054883

  * [6] L. Amendola, et al., Cosmology and fundamental physics with the Euclid satellite, Living Rev. Rel. 21 (1) (2018) 2.  arXiv:1606.00180, doi:10.1007/s41114-017-0010-3.
  * [7] H. E. S. Velten, R. F. vom Marttens, W. Zimdahl, Aspects of the cosmological “coincidence problem”, Eur. Phys. J. C74 (11) (2014) 3160.  arXiv:1410.2509, doi:10.1140/epjc/s10052-014-3160-4.
  * [8] I. Zlatev, L. Wang, P. J. Steinhardt, Quintessence, cosmic coincidence, and the cosmological constant, Phys. Rev. Lett. 82 (1999) 896–899.  arXiv:astro-ph/9807002, doi:10.1103/PhysRevLett.82.896.
  * [9] S. Weinberg, The cosmological constant problem, Rev. Mod. Phys. 61 (1989) 1–23.  doi:10.1103/RevModPhys.61.1.
  * [10] S. Weinberg, Theories of the cosmological constant, in: Critical dialogues in cosmology, 1996, pp. 195–203.  arXiv:astro-ph/9610044.
  * [11] Observation of gravitational waves from a binary black hole merger, Phys. Rev. Lett. 116 (2016) 061102.  doi:10.1103/PhysRevLett.116.061102.
URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102

  * [12] C. de Rham, S. Melville, Gravitational rainbows: Ligo and dark energy at its cutoff, Phys. Rev. Lett. 121 (2018) 221101\.  doi:10.1103/PhysRevLett.121.221101.
URL https://link.aps.org/doi/10.1103/PhysRevLett.121.221101

  * [13] P. Creminelli, F. Vernizzi, Dark energy after gw170817 and grb170817a, Phys. Rev. Lett. 119 (2017) 251302.  doi:10.1103/PhysRevLett.119.251302.
URL https://link.aps.org/doi/10.1103/PhysRevLett.119.251302

  * [14] J. Lesgourgues, S. Pastor, Massive neutrinos and cosmology, Phys. Rept. 429 (2006) 307–379.  arXiv:astro-ph/0603494, doi:10.1016/j.physrep.2006.04.001.
  * [15] M. G. Betti, et al., Neutrino physics with the PTOLEMY project: active neutrino properties and the light sterile case, JCAP 1907 (2019) 047.  arXiv:1902.05508, doi:10.1088/1475-7516/2019/07/047.
  * [16] M. Takada, E. Komatsu, T. Futamase, Cosmology with high-redshift galaxy survey: neutrino mass and inflation, Phys. Rev. D73 (2006) 083520.  arXiv:astro-ph/0512374, doi:10.1103/PhysRevD.73.083520.
  * [17] R. Jimenez, C. P. Garay, L. Verde, Neutrino footprint in Large Scale Structure, Phys. Dark Univ. 15 (2017) 31–34.  arXiv:1602.08430, doi:10.1016/j.dark.2016.11.004.
  * [18] M. Moresco, R. Jimenez, L. Verde, A. Cimatti, L. Pozzetti, C. Maraston, D. Thomas, Constraining the time evolution of dark energy, curvature and neutrino properties with cosmic chronometers, JCAP 1612 (12) (2016) 039.  arXiv:1604.00183, doi:10.1088/1475-7516/2016/12/039.
  * [19] C. Wetterich, Growing neutrinos and cosmological selection, Physics Letters B 655 (5-6) (2007) 201–208.  arXiv:0706.4427, doi:10.1016/j.physletb.2007.08.060.
  * [20] C. Y. Cardall, G. M. Fuller, Neutrino oscillations in curved spacetime: A heuristic treatment, Phys. Rev. D 55 (1997) 7960–7966.  doi:10.1103/PhysRevD.55.7960.
URL https://link.aps.org/doi/10.1103/PhysRevD.55.7960

  * [21] L. Buoninfante, G. G. Luciano, L. Petruzziello, L. Smaldone, Neutrino oscillations in extended theories of gravity, Phys. Rev. D 101 (2) (2020) 024016\.  arXiv:1906.03131, doi:10.1103/PhysRevD.101.024016.
  * [22] M. Blasone, G. Lambiase, G. G. Luciano, L. Petruzziello, L. Smaldone, Time-energy uncertainty relation for neutrino oscillations in curved spacetime, Class. Quant. Grav. 37 (15) (2020) 155004.  arXiv:1904.05261, doi:10.1088/1361-6382/ab995c.
  * [23] M. Blasone, P. Jizba, L. Smaldone, Flavor Energy uncertainty relations for neutrino oscillations in quantum field theory, Phys. Rev. D 99 (1) (2019) 016014\.  arXiv:1810.01648, doi:10.1103/PhysRevD.99.016014.
  * [24] A. W. Brookfield, C. van de Bruck, D. F. Mota, D. Tocchini-Valentini, Cosmology of mass-varying neutrinos driven by quintessence: Theory and observations, Phys. Rev. D 73 (2006) 083515.  doi:10.1103/PhysRevD.73.083515.
URL https://link.aps.org/doi/10.1103/PhysRevD.73.083515

  * [25] R. Fardon, A. E. Nelson, N. Weiner, Dark energy from mass varying neutrinos, Journal of Cosmology and Astroparticle Physics 2004 (10) (2004) 005–005.  doi:10.1088/1475-7516/2004/10/005.
URL https://doi.org/10.1088%2F1475-7516%2F2004%2F10%2F005

  * [26] D. B. Kaplan, A. E. Nelson, N. Weiner, Neutrino oscillations as a probe of dark energy, Phys. Rev. Lett. 93 (2004) 091801.  doi:10.1103/PhysRevLett.93.091801.
URL https://link.aps.org/doi/10.1103/PhysRevLett.93.091801

  * [27] H. Mohseni Sadjadi, H. Yazdani Ahmadabadi, Damped neutrino oscillations in a conformal coupling model, Phys. Rev. D 103 (2021) 065012.  doi:10.1103/PhysRevD.103.065012.
URL https://link.aps.org/doi/10.1103/PhysRevD.103.065012

  * [28] H. Mohseni Sadjadi, A. P. Khosravi, Symmetry breaking, and the effect of matter density on neutrino oscillation, JCAP 04 (2018) 008.  arXiv:1711.06607, doi:10.1088/1475-7516/2018/04/008.
  * [29] A. Esmaili, Y. Farzan, Implications of the Pseudo-Dirac Scenario for Ultra High Energy Neutrinos from GRBs, JCAP 12 (2012) 014.  arXiv:1208.6012, doi:10.1088/1475-7516/2012/12/014.
  * [30] J. F. Beacom, N. F. Bell, D. Hooper, J. G. Learned, S. Pakvasa, T. J. Weiler, Pseudo-dirac neutrinos: A challenge for neutrino telescopes, Phys. Rev. Lett. 92 (2004) 011101\.  doi:10.1103/PhysRevLett.92.011101.
URL https://link.aps.org/doi/10.1103/PhysRevLett.92.011101

  * [31] S. Ando, M. Kamionkowski, I. Mocioiu, Neutrino oscillations, $\mathrm{\text{lorentz}}/cpt$ violation, and dark energy, Phys. Rev. D 80 (2009) 123522.  doi:10.1103/PhysRevD.80.123522.
URL https://link.aps.org/doi/10.1103/PhysRevD.80.123522

  * [32] N. Klop, S. Ando, Effects of a neutrino-dark energy coupling on oscillations of high-energy neutrinos, Phys. Rev. D 97 (2018) 063006.  doi:10.1103/PhysRevD.97.063006.
URL https://link.aps.org/doi/10.1103/PhysRevD.97.063006

  * [33] F. Simpson, R. Jimenez, C. Pena-Garay, L. Verde, Dark energy from the motions of neutrinos, Physics of the Dark Universe 20 (2018) 72.  arXiv:1607.02515, doi:10.1016/j.dark.2018.04.002.
  * [34] A. R. Khalifeh, R. Jimenez, Spinors and Scalars in curved spacetime: Neutrino dark energy (DEν), Phys. Dark Univ. 31 (2021) 100777.  arXiv:2010.08181, doi:10.1016/j.dark.2021.100777.
  * [35] G. W. Horndeski, Second-order scalar-tensor field equations in a four-dimensional space, Int. J. Theor. Phys. 10 (1974) 363–384.  doi:10.1007/BF01807638.
  * [36] T. Kobayashi, Horndeski theory and beyond: a review, Rept. Prog. Phys. 82 (8) (2019) 086901.  arXiv:1901.07183, doi:10.1088/1361-6633/ab2429.
  * [37] B. Ratra, P. J. E. Peebles, Cosmological consequences of a rolling homogeneous scalar field, Phys. Rev. D 37 (1988) 3406–3427.  doi:10.1103/PhysRevD.37.3406.
URL https://link.aps.org/doi/10.1103/PhysRevD.37.3406

  * [38] S. Tsujikawa, Quintessence: A Review, Class. Quant. Grav. 30 (2013) 214003.  arXiv:1304.1961, doi:10.1088/0264-9381/30/21/214003.
  * [39] N. Birrell, P. Davies, Quantum Fields in Curved Space, chapter 3.8, Cambridge Monographs on Mathematical Physics, Cambridge Univ. Press, Cambridge, UK, 1984\.  doi:10.1017/CBO9780511622632.
  * [40] M. Lanzagorta, Quantum Information in Gravitational Fields, chapter 5, 2053-2571, Morgan & Claypool Publishers, 2014.  doi:10.1088/978-1-627-05330-3.
URL http://dx.doi.org/10.1088/978-1-627-05330-3

  * [41] V. Mukhanov, S. Winitzki, Introduction to quantum effects in gravity, Cambridge University Press, 2007.
  * [42] A. Capolupo, G. Lambiase, A. Quaranta, Neutrinos in curved spacetime: Particle mixing and flavor oscillations, Phys. Rev. D 101 (2020) 095022.  doi:10.1103/PhysRevD.101.095022.
URL https://link.aps.org/doi/10.1103/PhysRevD.101.095022

  * [43] S. M. Carroll, Spacetime and Geometry, Cambridge University Press, 2019.
  * [44] S. P. Mikheyev, A. Y. Smirnov, Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos, Sov. J. Nucl. Phys. 42 (1985) 913–917.
  * [45] S. P. Mikheev, A. I. Smirnov, Resonant amplification of $\nu$ oscillations in matter and solar-neutrino spectroscopy., Nuovo Cimento C Geophysics Space Physics C 9 (1986) 17–26.  doi:10.1007/BF02508049.
  * [46] C. Giunti, C. W. Kim, Fundamentals of Neutrino Physics and Astrophysics, chapter9.3, Oxford University Press, 2007.
  * [47] S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek, A. S. Bolton, J. R. Brownstein, A. Burden, C.-H. Chuang, J. Comparat, A. J. Cuesta, K. S. Dawson, D. J. Eisenstein, S. Escoffier, H. Gil-Marín, J. N. Grieb, N. Hand, S. Ho, K. Kinemuchi, D. Kirkby, F. Kitaura, E. Malanushenko, V. Malanushenko, C. Maraston, C. K. McBride, R. C. Nichol, M. D. Olmstead, D. Oravetz, N. Padmanabhan, N. Palanque-Delabrouille, K. Pan, M. Pellejero-Ibanez, W. J. Percival, P. Petitjean, F. Prada, A. M. Price-Whelan, B. A. Reid, S. A. Rodríguez-Torres, N. A. Roe, A. J. Ross, N. P. Ross, G. Rossi, J. A. Rubiño-Martín, S. Saito, S. Salazar-Albornoz, L. Samushia, A. G. Sánchez, S. Satpathy, D. J. Schlegel, D. P. Schneider, C. G. Scóccola, H.-J. Seo, E. S. Sheldon, A. Simmons, A. Slosar, M. A. Strauss, M. E. C. Swanson, D. Thomas, J. L. Tinker, R. Tojeiro, M. V. Magaña, J. A. Vazquez, L. Verde, D. A. Wake, Y. Wang, D. H. Weinberg, M. White, W. M. Wood-Vasey, C. Yèche, I. Zehavi, Z. Zhai, G.-B. Zhao, The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Monthly Notices of the Royal Astronomical Society 470 (3) (2017) 2617–2652.  arXiv:https://academic.oup.com/mnras/article-pdf/470/3/2617/18315003/stx721.pdf, doi:10.1093/mnras/stx721.
URL https://doi.org/10.1093/mnras/stx721

  * [48] S. Dodelson, F. Schmidt, Modern Cosmology, 2nd edition, Academic Press, 2020.
  * [49] Review of particle physics, Phys. Rev. D 98 (2018) 030001.  doi:10.1103/PhysRevD.98.030001.
URL https://link.aps.org/doi/10.1103/PhysRevD.98.030001

  * [50] The detection of a sn iin in optical follow-up observations of icecube neutrino events, The Astrophysical Journal 811 (1) (2015) 52.  doi:10.1088/0004-637x/811/1/52.
URL https://doi.org/10.1088/0004-637x/811/1/52

  * [51] C. Wetterich, Cosmology and the Fate of Dilatation Symmetry, Nucl. Phys. B 302 (1988) 668–696.  arXiv:1711.03844, doi:10.1016/0550-3213(88)90193-9.
  * [52] Y. Fujii, Origin of the gravitational constant and particle masses in a scale-invariant scalar-tensor theory, Phys. Rev. D 26 (1982) 2580–2588.  doi:10.1103/PhysRevD.26.2580.
URL https://link.aps.org/doi/10.1103/PhysRevD.26.2580

  * [53] L. H. Ford, Cosmological-constant damping by unstable scalar fields, Phys. Rev. D 35 (1987) 2339–2344.  doi:10.1103/PhysRevD.35.2339.
URL https://link.aps.org/doi/10.1103/PhysRevD.35.2339

  * [54] M. Fukugita, S. Watamura, M. Yoshimura, Light pseudoscalar particle and stellar energy loss, Phys. Rev. Lett. 48 (1982) 1522–1525.  doi:10.1103/PhysRevLett.48.1522.
URL https://link.aps.org/doi/10.1103/PhysRevLett.48.1522

  * [55] C. P. Burgess, J. M. Cline, New class of majoron-emitting double-$\beta$ decays, Phys. Rev. D 49 (1994) 5925–5944.  arXiv:hep-ph/9307316, doi:10.1103/PhysRevD.49.5925.
  * [56] E. Witten, Lepton number and neutrino masses, Nucl. Phys. B Proc. Suppl. 91 (2001) 3–8.  arXiv:hep-ph/0006332, doi:10.1016/S0920-5632(00)00916-6.
  * [57] Y. Chikashige, R. N. Mohapatra, R. D. Peccei, Are There Real Goldstone Bosons Associated with Broken Lepton Number?, Phys. Lett. B 98 (1981) 265–268.  doi:10.1016/0370-2693(81)90011-3.
  * [58] D. Chowdhury, J. Martin, C. Ringeval, V. Vennin, Assessing the scientific status of inflation after Planck, Phys. Rev. D 100 (8) (2019) 083537.  arXiv:1902.03951, doi:10.1103/PhysRevD.100.083537.
  * [59] A. Gando, Y. Gando, H. Hanakago, H. Ikeda, K. Inoue, R. Kato, M. Koga, S. Matsuda, T. Mitsui, T. Nakada, et al., Limits on majoron-emitting double-$\beta$decays of136xe in the kamland-zen experiment, Physical Review C 86 (2) (Aug 2012).  doi:10.1103/physrevc.86.021601.
URL http://dx.doi.org/10.1103/PhysRevC.86.021601

  * [60] K. J. Kelly, P. A. Machado, Multimessenger astronomy and new neutrino physics, Journal of Cosmology and Astroparticle Physics 2018 (10) (2018) 048–048.  doi:10.1088/1475-7516/2018/10/048.
URL http://dx.doi.org/10.1088/1475-7516/2018/10/048

  * [61] H. Athar, J. F. Nieves, Matter effects on neutrino oscillations in gravitational and magnetic fields, Phys. Rev. D 61 (2000) 103001.  arXiv:hep-ph/0001069, doi:10.1103/PhysRevD.61.103001.
  * [62] M. Lei, N. Steinberg, J. D. Wells, Probing Non-Standard Neutrino Interactions with Supernova Neutrinos at Hyper-K, JHEP 01 (2020) 179.  arXiv:1907.01059, doi:10.1007/JHEP01(2020)179.
  * [63] M. G. Aartsen, et al., IceCube-Gen2: The Window to the Extreme Universe, J. Phys. G 48 (6) (2021) 060501.  arXiv:2008.04323, doi:10.1088/1361-6471/abbd48.
  * [64] S. Böser, M. Kowalski, L. Schulte, N. L. Strotjohann, M. Voge, Detecting extra-galactic supernova neutrinos in the Antarctic ice, Astropart. Phys. 62 (2015) 54–65.  arXiv:1304.2553, doi:10.1016/j.astropartphys.2014.07.010.
  * [65] S. Capozziello, M. De Laurentis, Extended Theories of Gravity, Phys. Rept. 509 (2011) 167–321.  arXiv:1108.6266, doi:10.1016/j.physrep.2011.09.003.
  * [66] R. Maartens, K. Koyama, Brane-World Gravity, Living Rev. Rel. 13 (2010) 5.  arXiv:1004.3962, doi:10.12942/lrr-2010-5.