Distinguishing Dark Energy Models with Neutrino Oscillations

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\^itre-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 in future underground, terrestrial or extraterrestrial neutrino telescopes, which would assess probing Dark Energy models with this technique.


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/CP T violating theories [31,32], where CP T 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 Λ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 = c = 1, where c is the speed of light and 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.

Generalized Formalism
Let us consider a general interaction of a spinor field, ψ, with a classical scalar field, ϕ, in curved spacetime with metric g µν . This is described by the following action: In here, g is the determinant of g µν and R is the Ricci scalar, the trace of the Ricci tensor R µν . Moreover, D µ 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, D µ reduces to ∂ µ , 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 (ϕ), which describes self interactions of ϕ. Furthermore, γ µ = {γ 0 , γ 1 , γ 2 , γ 3 } are the four Dirac matrices,ψ = ψ † γ 0 , with ψ † being the complex transpose of ψ, and m is the spinor field's mass. Finally, ζ is the coupling constant for the general interaction term between ψ and ϕ, Θ ψ,ψ, ϕ, X µ ψ , X μ ψ , X µ ϕ , with X µ ψ X μ ψ = D µ ψ D µψ and X µ ϕ = ∂ µ ϕ. By setting the variational derivative of the action (1) with respect to(w.r.t) ψ to 0, we get the modified Dirac equation for ψ, i.e.
At this stage, one could say that the interaction is of the most general form, for neither Θ 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 Θ w.r.tψ. 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: 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 ξ. On the other hand, the second term is a flavor-specific term, with the coupling strength ξ f depending on which flavor is being considered. Another term that could be added is a kinetic coupling, such as ξF (ϕ, X µ ϕ )γ µ D µ ψ. 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 ΛCDM and quintessence(in the context of DE ν ), it turns out that having is more useful, and will therefore be used in the following sections

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.

Transition Amplitude's Evolution
The first step in studying neutrino oscillations is to expand a state of flavor α, |ν α , in terms of mass eigenstates, |ν j : where λ is the monotonically increasing affine parameter along the neutrino world-line and U αj is the two-flavor mixing matrix, given by: with θ being the mixing angle. Moreover, Φ(λ) is the mass eigenstate's evolution operator [20]: where P µ is the 4-momentum operator, dx µ /dλ is a null vector tangent to the neutrino world-line, and λ 0 (λ) 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 and therefore the transition amplitude between states α and β, The ultimate goal is to find the transition probability between flavors α and β, i.e.
For the purpose at hand, let us start with a system of two-flavors, electron and muon neutrinos ν e and ν µ , respectively, that is α, β = e, µ. Let be a vector of spinor fields. The modified Dirac equation (2) for this system becomes: where the vacuum mass matrix in flavor space is given by with m 1 , m 2 being the masses of mass states |ν 1 , |ν 2 , respectively. Now we introduce the explicit form of the covariant derivative D µ . 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 γ a , and link it to the general one using tetrad fields e µ a , where Greek indices correspond to general coordinates, while Latin ones for the local system. With this, we can write where is called the spin-connection which takes into account the gravitational effect on the particle's spin. In eq. (16), [γ b , γ c ] is the commutator of γ a and γ b , and ∇ µ is the usual covariant derivative of General Relativity [43]. With these relations, one can then show that with where 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: 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: where It should be noted here thatm is not a massvarying 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 which finally implies, from eqs. (7) and (10), that with V I = −A µ dx µ /dλ. 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λ, and second, P i and dx i /dλ are assumed parallel, that is

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 where I is the 2 × 2 identity matrix in flavor space and∆ 2 . 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 ν e state, for instance, the evolution equation for the transition amplitudes becomes: Notice that the gravitational contribution A Gµ 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 M , which has (26) as eigenvalues, andŨ = cosθ sinθ − sinθ cosθ (27) as the unitary matrix that diagonalizes it, with 2 In analogy with the flavor-mass bases transformation, let us define as a vector of transition amplitudes in an effective mass basis, {ν − , ν + }, that takes into account the neutrino interaction with gravity and DE. Using the unitarity ofŨ and eq. (25), it can be shown that φ e satisfies: 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θ/dλ = 0, and the transition amplitudes evolve as: for j = +, −, where φ ej (0) is the initial condition and 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θ from eq. (28) w.r.t λ: Since we expect gravitational and DE effects to be small compared to the vacuum oscillations, that is ξ, ξ f 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λ = 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 ν e state, we can write where, by construction, an initial ν e state corresponds to Ψ ee (0) = 1 and Ψ eµ (0) = 0. By acting with the inverse transformation of eq. (29), we can calculate the amplitude Ψ eµ , and thus, with the initial conditions eq. (34), we finally obtain the ν e → ν µ transition probability: Let us now look in more detail into the oscillating term in eq. (35). If we use the above mentioned approximation (ξ, ξ f 1), one can show that 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λ = dt/E = dx/E, we get where L is the distance traveled by the neutrinos. This results in 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θ 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θ from eq. (28) into eq. (35), and then expand all functions of the interactions up to 1 st order, we get where P vac νe→νµ = sin 2 2θ sin 2 ω vac is the transition probability in vacuum, with frequency and 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.

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 Λ and Quintessence.

Λ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, ΛCDM has two main components in the late universe: a cosmological constant DE, Λ, and cold Dark Matter(CDM). The metric of spacetime that best describes it is FLRW: where t is cosmic time, x i , for i = 1, 2, 3, are comoving spatial coordinates, δ 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: where H =ȧ/a is the Hubble parameter, with H 0 being its value today and Ω m0 , Ω Λ0 are today's matter and Λ 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 Λ [48]. To see the effect of DE in this model on neutrino oscillations, let us start by noting that when DE is Λ, ξ = ξ f = 0 in eq. (4), implying the automatic satisfaction of the adiabaticity condition, dθ/dλ = 0, from eq. (33), and thus sinθ(cosθ) = sin θ(cos θ). Therefore, from the first equality in eq. (36), we define with the second equality meaning eq. (36) when DE is Λ. Here, t 0 is today, t em is the time of neutrino emission and E = dt/dλ is the 0 th component of the null tangent vector dx µ /dλ, 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), where a em (z em ) is the scale factor(redshift) at neutrino emission, and with the usual normalization a 0 = 1and 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, and E = E 0 (1 + z e ), respectively, we get, Finally, inserting eqs. (46) and (48) in eq. (35) gives the two-flavor neutrino oscillation probability in the ΛCDM model, and in the SA, P SA = sin 2 2θ sin 2 ω SA .
In Figure 1, we plot the evolution of P Λ (solid black curve) and P SA (dotted blue curve) as a function of redshift, to compare the two approaches. To this end, we took ∆ 2 m = 7.53×10 −5 eV 23 and E 0 = 10 16 eV, a value to which neutrino detectors are on average sensitive to [50]. Further, we used Ω m0 = 0.315, Ω Λ0 = 0.685 and H 0 = 1.44 × 10 −33 eV as reported in [3]. For redshifts higher than ∼ 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 ≤ z ≤ 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 3 Here we used mass states 1 and 2 from Ref. [49] as ours. One can check that other values of ∆ 2 m 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 −1 0 ∼ O(10 33 )eV, which wipes out the O(10 2 )eV 2 difference between ∆ 2 m s. 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 Λ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.

Quintessence
As a homogeneous canonical scalar field minimally coupled to gravity, Quintessence 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 Λ 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 µ = ∂ µ ϕ 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 whereφ = ∂ t ϕ and V ϕ (ϕ) is the potential energy of ϕ. Therefore, eq. (44) becomes where 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: For ϕ to produce an accelerated expansion, it should satisfy the condition: where ρ Λ0 is the energy density of a cosmological constant today. This means that, first, in eq. (52), Ω ϕ ≈ Ω Λ0 , and second, we can write 4 4 The 7/2 factor is to reduce numerical factors clustering.
where < ρ Λ0 ∼ O(10 −11 eV 4 ) 5 , and thus, from eq. (54), we geṫ On the neutrino's side, this type of interaction results in from which one can show, using eqs. (26), (28) and (33), that and respectively. To check if the adiabaticity condition is satisfied for the current case, differentiate eq. (58) w.r.t λ and insert it in eq. (61), to find From the fact that H 0 ∼ O(10 −33 )eV [48], E 0 ∼ O(10 16 )eV(typical value for high-energy neutrinos [50]), ξ e,µ ∼ O(10 −14 )eV −1 [34] and ∼ O(10 −11 )eV 4 , one can see that dθ/dλ v ± , and therefore the adiabaticity condition still holds, resulting in an oscillation frequency Finally, from eq. (35), the two-flavor oscillation probability in the case when DE is quintessence is: To study the difference between this model and ΛCDM, we plot(figure 3) eq.(64) for values of ξ i , i = e, µ, ranging from 10 −17 to 10 −14 eV −1 , in addition to eq.(49) for the ΛCDM case. We have checked that smaller values of ξ i do not produce any noticeable deviation from ΛCDM, while already at 10 −14 we can see from figure 3 that the deviation is ∼ 50% at z = 2. That is the reason why we focus on this range of values of the couplings ξ i . Moreover, we use the same parameters used to produce figures 1 and 2(see text after eq. (50)), in addition to cos 2θ = 0.4 [49]. As the strength of the coupling increases, the difference between the two models starts to become apparent at redshift ∼ 0.5, which is expected since then DE is becoming more dynamical than in the case of Λ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 and ξ f of this particular DE model. First, if quintessence is slow-rolling, but not ultra slow-rolling, then ε(see footnote 5) cannot be too small [58]. Taking ∈ [10 −14 , 10 −12 ], which corresponds to ε ∈ [10 −3 , 10 −1 ], we find that ξ f ∈ [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 ∼ O(10 −13 ), for instance, we find that the difference between quintessence and ΛCDM starts to become appreciable(i.e. more that a few %) for ξ f ∼ O(10 −14 − 10 −13 ). This is still consistent with Particle Physics constraints for this type of coupling, which is ξ f 10 −7 eV −1 [59,60]. Moreover, for values of ξ 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 ΛCDM, the neutrino-quintessence interaction is affected by the value of ∆ 2 m 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 ∆ 2 m = 2.51 × 10 −3 eV 2 (normal hierarchy), while in the lower panel we use ∆ 2 m = −2.56 × 10 −3 eV 2 (inverted hierarchy), with cos 2θ ∼ 0.2 for both. These values of ∆ 2 m 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 Λ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(∼ 10 −14 ), as is apparent from the two panels of figure 4. Therefore, such a neutrinoquintessence interaction could require some fine-tuning to match future obser- vations, which puts it at equal, or less, footing with ΛCDM 6 .

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 νe can be expressed as [61,62]: where F 0 να is the flux of neutrinos with flavor α 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 να→νe . More specifically to our case, depending on which DE model is considered, eqs. (49) and (64) will give different P να→ν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∼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].

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 Λ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, Λ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 ν e from a type Ia supernova(SN), for example, and we want to calculate their flux (which depends on the ν e 's survival probability), SA would give a result ∼ 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, Λ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 ∆ 2 m 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 Λ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 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.

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).