{'fraction_non_alphanumeric': 0.03391253684282679, 'fraction_numerical': 0.011464116800329015, 'mean_word_length': 4.3959315765141005, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Plane-wave and wavepacket neutrino flavor oscillations in vacuum in conformal coupling models

We investigate neutrino flavor oscillations in vacuum within a general conformal coupling model. We first examine the flavor oscillations within the plane-wave description of neutrinos, and then we extend our analysis to the wavepacket-based formalism. In both cases, we derive the general formulas for the flavor transition probability in arbitrary static and spherically symmetric spacetimes. We thoroughly discuss and assess in both cases the different possible ways -- dictated by the presence of the conformal coupling -- of computing the flavor transition probability. We show that the conformal invariance of the Dirac equation implies that the effect of conformal coupling on neutrino flavor oscillations modifies the oscillation length, but preserves the coherence length and unitarity. A detailed application to two-flavor neutrinos is then made and a numerical analysis is conducted within the well-known chameleon conformal coupling model.

We investigate neutrino flavor oscillations in vacuum within a general conformal coupling model. We first examine the flavor oscillations within the plane-wave description of neutrinos, and then we extend our analysis to the wavepacket-based formalism. In both cases, we derive the general formulas for the flavor transition probability in arbitrary static and spherically symmetric spacetimes. We thoroughly discuss and assess in both cases the different possible ways -dictated by the presence of the conformal coupling -of computing the flavor transition probability. We show that the conformal invariance of the Dirac equation implies that the effect of conformal coupling on neutrino flavor oscillations modifies the oscillation length, but preserves the coherence length and unitarity. A detailed application to two-flavor neutrinos is then made and a numerical analysis is conducted within the well-known chameleon conformal coupling model.

I. INTRODUCTION
Attempts to explain the dark matter and dark energy problems of modern cosmology [1] have led to various proposals for modifying the Einstein-Hilbert gravitational action. Such proposals can be classified into two main categories. In the first, one tweaks the gravitational action by putting in extra invariant curvature terms besides the Ricci scalar R without adding any extra fields. In the second, one inserts into the action extra fields and geometric degrees of freedom in addition to possible extra invariant curvature terms. The latter models are known as scalar-tensor theories of gravity. Among the various proposals in the latter category are the conformal coupling models in which an extra field added to the action is conformally coupled to the matter sector of the action. The models that attracted much attention are the so-called chameleon model [2] and the symmetron model [3] as they exhibit a screening of the unwanted fifth force induced by the scalar field on matter. Although screening is also realized via other mechanisms, such as the Vainshtein mechanism [4,5], and within dilaton models [6], we focus in this paper on the chameleon-and symmetron-type models as they are the most studied in the literature [7].
The gravitational sector in those conformal coupling models contains only the Einstein-Hilbert action and the usual action of a scalar field φ(x) with a kinetic term and a potential term. However, the matter sector of the models displays a conformal coupling between the gravitational field and the matter fields Ψ i (x). Such a gravitational coupling is said to be conformal because it amounts to replacing the curved-spacetime matter Lagrangian L m (Ψ i , g µν ) by the Lagrangian L m (Ψ i ,g µν ).
The metricg µν is related to the spacetime metric g µν of the gravitational sector by a Weyl conformal transformation [8] g µν = Ω 2 (φ)g µν , where the nonvanishing conformal factor Ω(φ) is a functional of the scalar field φ(x). In the chameleon model the conformal factor has the form Ω(φ) = exp(β φ/M p ), where M p stands for the reduced Planck mass and β is an arbitrary dimensionless constant [2]. In the symmetron model the conformal factor has the form Ω(φ) = 1 + 1 2 φ 2 /M 2 [3], where M is an arbitrary mass scale. In both models, one assumes, for simplicity, a universal conformal coupling in the sense that the scalar field φ couples identically to all matter fields. For more generality, however, one may introduce a different functional Ω i (φ) and/or a different scalar field φ i for each matter field Ψ i . Indeed, this would offer one extra possibility for violating the Equivalence Principle [9].
On the other hand, the well-known theoretical proposal for neutrino flavor transition [10,11] and its experimental discovery [12][13][14][15][16][17][18][19][20][21][22][23] have given rise to a wide variety of scenarios for the origin (and even possible variation) of neutrino masses [24]. As the masses of neutrinos are the ones responsible for giving rise to the observed neutrino flavor oscillations, any possible variation of the masses could be confirmed by measuring neutrino oscillations [25][26][27][28]. However, among the various proposed mechanisms that are of interest to us here are those that lead to varying masses for neutrinos based on the interaction of the latter with a scalar field [29][30][31][32][33][34] and, more specifically, the conformal coupling scenario. The effect of this kind of possible variation of neutrinos' masses on neutrino oscillations has more recently been examined by various authors [35][36][37][38][39]. Now, it is also well known that the Dirac equation (and hence the fermion Lagrangian) in curved spacetime is conformally invariant [40]. Therefore, replacing the spacetime metric g µν by the conformal metricg µν inside the matter Lagrangian simply amounts to using the original metric g µν after replacing the fermion's mass m, its energy E and its wavefunction ψ by a new massm, a new energyẼ and a new wavefunctionψ. Thus, when a free neutrino mass eigenstate of mass m j and of energy E j is conformally coupled to the spacetime, it actually propagates in the curved spacetime of metric g µν = Ω −2g µν as a 'free' neutrino mass eigenstate of massm j and of energyẼ j . Under such a conformal deformation of the metric, the masses and energies of quantum particles are related, as usual [41], bym j = Ωm j andẼ j = ΩE j , respectively. Therefore, since the characteristic length L jk osc of neutrino flavor oscillations is 4πĒ/|m 2 j −m 2 k | (whereĒ is the average energy of the neutrinos and j and k stand for the two mass eigenstates 1 ), we expect that the conformal coupling would simply lead to a new oscillation lengthL jk osc given by 4πĒ/|m 2 j −m 2 k |. In other words, we expect that a conformal coupling of neutrinos would lead to a new oscillation length which is roughly of the formL jk osc ∼ L jk osc /Ω. Given that the oscillation phase without conformal coupling is proportional to L p /L osc [42], where L p is the proper distance traveled by the conformally coupled neutrinos, the expected effect of the conformal coupling is then to induce oscillation phases that are proportional to the conformal factor.
The previous reasoning applies to neutrinos viewed as plane waves of definite wavelengths. Nevertheless, neutrino flavor oscillations are also very well studied in the literature using the wavepacket formalism as required by the uncertainties associated with the neutrinos production and detection processes [42]. In Refs. [43,44], the effect of curved spacetime on the flavor oscillations within the wavepacket formalism has been investigated in much detail and shown to give rise to an oscillation phase that contains the ratio I/L jk osc , where I is an integral involving the components of the metric. In addition, however, the oscillation phase contains extra metricdependent terms, as well as a damping term that also depends in a nontrivial way on the metric components. We therefore expect that while the plane-wave formalism might lead to simple expressions for the effect of a conformal coupling on neutrino oscillations, the wavepacket formalism could potentially display a more involved expression for the effect. Furthermore, the wavepacket approach in curved spacetime brings with it an additional complication due to the two possible ways of dealing with the propagation of the mass eigenstates wavepackets [44].
Our first goal in this paper is to make the above handwaving reasoning more precise by systematically deriving the flavor transition probability for conformally coupled neutrinos within both the plane-wave approach and the wavepacket approach. Our second goal is to give a physical interpretation of the subtleties brought in by the conformal coupling when dealing with wavepackets.
The remainder of this paper is organized as follows. In Sec. II, we derive the flavor transition probability for conformally coupled neutrinos within the plane-wave formalism. In Sec. III, we briefly recall the wavepacket approach to neutrino oscillations and the various subtleties it conceals, and then apply the approach to conformally coupled neutrinos. In Sec. IV, we apply to the case of two-flavor neutrinos the general transition probabilities obtained in sections II and III. Finally, although we consider throughout this paper the case of universal conformal couplings, Sec. V is devoted to discussing the way the results of sections II and III generalize further to the case of non-universal conformal couplings. We conclude this pa-1 See, e.g., Ref. [42] for a rich introduction to neutrino oscillations and for many relevant references.
per with a brief section in which we summarise and discuss our various results.

II. CONFORMALLY COUPLED PLANE-WAVE NEUTRINOS
The curved-spacetime Dirac equation describing the dynamics of a spinor wavefunction ψ j associated to a neutrino mass eigenstate of mass m j reads (iγ µ ∇ µ − m j ) ψ j = 0 2 . This equation is obtained from the flat-spacetime Dirac equation by replacing the gamma matrices γ a by the curvedspacetime gamma matrices γ µ = e µ a γ a and the partial derivatives ∂ µ by the spin covariant derivatives The spin connection ω ab µ is related to the vierbeins e a µ and the Christoffel symbols Γ λ µν by ω ab µ = e a ν ∂ µ e νb + Γ λ µν e a λ e νb . This curved-spacetime Dirac equation can be obtained from the fermion-matter Lagrangian [45] d where the metric determinant √ −g is that of the metric g µν , andψ is the adjoint spinor:ψ = ψ † γ 0 . However, for neutrinos conformally coupled to the metric via a scalar field φ their Lagrangian and their derived Dirac equation should both be written using the metricg µν = Ω 2 (φ)g µν , the curved-spacetime gamma matrices γ µ = Ω(φ)γ µ and the covariant derivatives∇ µ associated to the metricg µν . We shall assume φ and g µν are time-independent. As both the Lagrangian and the Dirac equation are conformally invariant, the conformally-coupled neutrinos' Lagrangian and their Dirac equation can be written in terms of the original spacetime metric g µν = Ω −2 (φ)g µν , respectively, as follows: Since in the new expressions of both the Lagrangian and the Dirac equation the metric used is the original g µν of the gravitational sector, the scalar field φ is now hiding inside the new effective massm j and the new wavefunctionψ j . Indeed, these are given in terms of the original mass m j and the original wavefunction ψ j , respectively, by [35,40]: The new effective massm j of each mass eigenstate becomes thus φ-dependent and, hence, a position-dependent mass; i.e., a varying mass. Note that the conformal transformation of ψ as given by Eq. (3) is consistent with the requirement of preserving unitarity. Indeed, when working with the metricg µν the probability density to be integrated over a spatial region of volumeṼ of space is √h ψ † ψ, whereh is the determinant of the induced spatial metric in the region. On the other hand, when working with the metric g µν , the probability density to be integrated over the spatial region of volume V is the product √ hψ †ψ which, according to the second identity in Eq. (3) and the fact that In other words, the probability density is conserved by conformally transforming both the spacetime and the Dirac spinor.
Within the metric g µν the conformally coupled mass eigenstates satisfy the mass-shell condition g µνp µ jp ν j = −m 2 j , where the four-momentump µ j carries the effect of the scalar field φ through the following identity: with p µ j = m j dx µ /dτ being the four-momentum that a neutrino mass eigenstate |ν j of mass m j would have in the spacetime of metric g µν if it were not conformally coupled.
Recall now that a neutrino mass eigenstate |ν j not conformally coupled to spacetime, emitted at a point (0, x A ) in a spacetime of metric g µν , is described at a detection point (t, x B ) of the spacetime by the following plane wave [42]: where the orthonormal basis states |ν j of the Hilbert space satisfy ν j |ν k = δ jk , and where the quantum phase Φ PW j is evaluated using the Stodolsky prescription for computing the accumulated quantum phases of freely propagating quantum particles in curved spacetimes [46]: Note that we have attached here the superscript "0" to the mass eigenstate |ν 0 j (t, x B ) and to its wavefunction ψ 0 j (t, x B ) in order to distinguish them from the corresponding quantities |ν j (t, x B ) and ψ j (t, x B ) of the conformally coupled neutrinos propagating within the metricg µν . Note also that we attached the superscript "PW" to the plane wave's quantum phase in order to distinguish it from the wavepacket's quantum phase Φ WP j of Sec. III. Therefore, according to our discussion in the previous paragraph and according to Eq. (4), a neutrino mass eigenstate |ν j conformally coupled to the scalar field φ within a spacetime of metric g µν can be described at the point (t, x B ) by a plane wave |ν j (t, x B ) that reads where the quantum phaseΦ PW j is given bỹ Note that the orthonormal basis states of the Hilbert space in Eq. (7) satisfy again ν j |ν k = δ jk . Before using the results (7) and (8) to discuss the flavor transitions, we shall pause here to discuss the following important observation. Our equation (7) differs from Eq. (30) of Ref. [35] mainly by the absence of any multiplicative factor made of powers of Ω(φ). Such an absence in Eq. (7) might seem to be at odds with the second identity in Eq. (3). One might indeed expect that the wavefunction in Eq. (7) would differ from the wavefunction in Eq. (5) only by the multiplicative factor Ω 3 2 (φ). One must bear in mind, however, that ψ j in Eq. (3) is the wavefunction of the conformally coupled neutrino mass eigenstate as it propagates within the metricg µν , whereas ψ 0 j in our equation (5) is the wavefunction of a neutrino mass eigenstate not conformally coupled to the spacetime. There is no such a link betweenψ j and ψ 0 j as the one given by Eq. (3) betweenψ j and ψ j . In fact, as we shall see shortly, the absence of any multiplicative factor made of powers of Ω(φ) in our equation (7) is what will make the total flavor transition probability add up to unity, in contrast to what one would have obtained if multiplicative factors appeared in Eq. (7) as in Ref. [35]. Now, a neutrino flavor |ν α is made out of a superposition of all three mass eigenstates |ν j via the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) unitary mixing matrix U * α j [42]: |ν α = j U * α j |ν j . Therefore, the quantum mechanical probability for a conformally coupled neutrino produced at point (0, x A ) as an α flavor to be detected at the point (t, x B ) as a β flavor is We have defined here, for convenience, the resulting phase difference by the usual notationΦ PW jk ≡Φ PW j −Φ PW k . It must be noted here that although the wavefunction involved in Eq. (9) is the conformally transformed one, the mixing matrix U * α j is not affected by the conformal coupling. The reason is that the elements of such a matrix are dimensionless constants that specify only how much of each mass eigenstate in the Hilbert space is contained in the superposition that gives rise to a given flavor state. The conformal coupling of the mass eigenstates affects only their individual masses and wavefunctions living in the configuration space, it does not affect the superposition of the basis states |ν j . A more formal presentation of this argument is given in Appendix A.
Eq. (9) is the general expression of the flavor transition probability of conformally coupled neutrinos propagating inside an arbitrary spacetime of metric g µν . In order to use this result, we shall consider for simplicity in the remainder of this paper a radial propagation of neutrinos inside a static and spherically symmetric spacetime. Let the spacetime metric of the conformal coupling model be static and spherically symmetric, of the form ds 2 = −A(r) dt 2 + B(r) dr 2 + r 2 dθ 2 + sin 2 θ dϕ 2 . (10) Here, the functions A(r) and B(r) are arbitrary functions of the radial coordinate r, which are everywhere regular except maybe at a singularity or on a horizon.
Note that although according to Birkoff's theorem one has A(r)B(r) = 1 for most physical cases of such spacetimes, we shall assume throughout this paper the more general case A(r)B(r) 1 in view of its application to neutrino oscillations within modified gravity theories [44,47] and in view of its application to neutrino oscillations within the gravitational field inside matter [44]. The usual case A(r)B(r) = 1 of the exterior Schwarzschild solution found within general relativity will be dealt with in Sec. IV.

A. Integrating independently over t and r
Within a spacetime of metric (10), the energyẼ j (p) and momentump j of a conformally coupled and radially propagating mass eigenstate of massm j are, as perceived by an observer at infinity, given by the following expressions: (11) In the first step of each equality we used Eq. (3), and in the second step we used Eq. (4). The spacetime (10) possesses a timelike Killing vector K ν = (1, 0, 0, 0), as well as a spacelike Killing vector R µ = (0, 0, 0, 1). The time-like Killing vector implies the conserved energies whereas the second Killing vector implies the conservation of the angular momentum which need not concern us here since we deal with radially propagating neutrinos. The conserved energy E j (p) is the energy an observer at infinity would perceive if the neutrino mass eigenstate of mass m j were freely propagating in the spacetime of metric g µν . Combining the identities in Eq. (11) with the mass-shell condition −m 2 On the other hand, plugging the second identity in Eq. (11) into the mass-shell condition (13), we arrive at In the second step, we expanded the expression up to the first order in the ratiom 2 j /Ẽ j (p). Plugging Eq. (14) into Eq. (8), and then integrating independently over t and r, the accumulated phase of each mass eigenstate at the detection point (t, x B ) takes the following form: In the second step, we used identities (3) and (12) that allowed us to move the conserved energy E j (p) and the constant mass m j out of the integrals, leaving behind the conformal factor Ω(φ) that gets thus integrated along the path of the neutrinos from the emission to the detection point. In the third step, we introduced, for convenience, the notations T , I and J for the three different integrals. The counterparts of these integrals in the case of neutrinos not conformally coupled are, respectively, the integrals T , I 1 and I 2 obtained in Ref. [44]. Inserting expression (15) of the phases into Eq. (9), we find the following T -dependent transition probability from an αflavor neutrino to a β-flavor neutrino: We have used here the notation E jk (p) to denote the difference E j (p) − E k (p), as well as the usual notation ∆m 2 jk to denote the difference m 2 j − m 2 k . We have now to consider the time-integral T to be unobservable and compute the average probability by averaging expression (16) over T . This option is what is usually adopted in the literature when computing the transition probability [42]. The reason for doing so is that while the distance traveled by the neutrinos might be known precisely, the times of emission and detection of the neutrinos might not. Averaging the transition probability (16) by integrating out the unobserved quantity T gives rise to the delta function δ[E j (p) − E k (p)], which thus simply eliminates the term iE jk (p)I inside the exponential and replaces E j (p) and E k (p) in the remaining two terms by the average energy E 0 the neutrinos would have if they were massless [42]. This leads then to the following expression for the averaged transition probability: where the algebraic oscillation length L jk osc is given by When Ω(φ) is a constant equal to unity, the general probability expression (17) reduces to the usual result for the case of neutrinos propagating in a curved spacetime without being conformally coupled. However, for the case of neutrinos conformally coupled and propagating within a Schwarzschild spacetime or within a Minkowski spacetime, for which A(r)B(r) = 1, our result (17) does not reduce to what is found in Ref. [35]. It is also worth stressing here the following important observation. The result (17) shows that the effect of the conformal coupling of neutrinos to the scalar field in screening models shows up as an interesting cumulative phenomenon that is able to take into account the variation of the scalar field with position. However, when the scalar field φ settles down to its minimum value φ 0 all along the neutrinos path, the factor Ω(φ) in integral J can be moved out of the integration symbol. This would then simply amount to keeping the same formulas for neutrino flavor oscillations in curved spacetimes in the absence of any conformal coupling, and to just redefine the squared masses m 2 j by the effective ones Ω(φ 0 )m 2 j . In other words, the fifth-force screening mechanism, for which such modified-gravity models were designed in the first place, also works for neutrino flavor oscillations. That is, within local dense environments such as our solar system, and as long as one relies on the neutrino masses inferred from local solarsystem neutrino-detection experiments, neutrino flavor oscillation in vacuum would not provide any signature of the conformal coupling in that case. Only outer-galactic neutrinos and neutrinos that traveled inside matter could possibly provide such a signature, for then the position dependence of the scalar field becomes significant. Now, the above procedure of averaging over time is amply justified in the absence of a conformal coupling of neutrinos. Indeed, the time-integral T just reduces then to the total unmeasured time T [44]. However, when neutrinos interact with a scalar field that displays a non-negligible variation with position along the path of the neutrinos towards the detector, actually performing the integral T becomes ambiguous. In fact, in such a case one can neither move the factor Ω(φ) out of the integral nor replace the factor Ω(φ) by its average value over the neutrinos path. In addition, as we shall see in Sec. III, the presence of the scalar field in the phase factor of neutrinos wavepackets makes the averaging-over-time approach even more problematic. For a rigorous treatment of the effect of the scalar field, one should therefore take into account the path of the neutrinos by using the link between dt and dr that follows by combining Eqs. (11) and (13). Showing how to properly implement such a procedure in detail is the subject of the next subsection.

B. Converting the time integral into a spatial integral
To take into account the presence of φ inside the time integral T when computing the phase, we need to express the time element dt in that integral in terms of the radial line element dr. Isolating dt/dτ in terms ofẼ j (p) from the first identity in Eq. (11) and then dividing the result by the dr/dτ that we isolate from Eq. (13), we find the following link between the coordinate elements dt and dr along the path of the j-th mass eigenstate: In the second step we expanded the square root up to the first order inm 2 j /Ẽ 2 j (p). We need now to use this link to find the quantum phase of each mass eigenstate. However, to achieve that purpose we cannot merely insert Eq. (19) into the time integral in Eq. (15) corresponding to the j-th mass eigenstate and then merely replace everywhere j by k for the case of the k-th mass eigenstate. The reason is that a neutrino is detected as having a certain flavor when the mass eigenstates making that flavor are all detected at the same point in space and at the same time. Therefore, since the detection position is r = r B for all the mass eigenstates, a common detection time is obtained only if one turns all time integrals into spatial integrals using a unique link between the elements dt and dr. The time integral that needs to be performed in Eq. (15) for any mass eigenstate should then be the time integral corresponding to the most massive state of the three mass eigenstates, for the wave front of the most massive mass eigenstate is the last one to arrive at r = r B .
Let us denote the mass and energy of the most massive state bym * andẼ * , respectively. Then the time integral that should be inserted into the calculation (15) of the phase for all the mass eigenstates is Inserting this time T * into the second line of Eq. (15), the resulting quantum mechanical phase of the j-th mass eigenstate is given by the following expressioñ Therefore, inserting into expression (9) of the transition probability the phase differenceΦ PW jk one obtains by subtracting from Eq. (21) the corresponding expression for the k-th mass eigenstate, the transition probability takes the following form Here, E jk (p) stands again for the difference E j (p) − E k (p). Note that there is no more time-dependence in this expression of the probability, and therefore no averaging over time needs to be performed within this approach.
Next, approximating the conserved energies E j (p) and E k (p) of each mass eigenstate by the common average energy E 0 of massless neutrinos, the transition probability (22) takes exactly the same form as the one obtained within the previous approach; i.e., it takes the form given by Eq. (17) in which the algebraic oscillation length L jk osc is again given by Eq. (18). Note that one may go up to the first order in m 2 /E and approximate the conserved energies E j (p) and E k (p) by writing E j (p) ≈ E 0 + ξm 2 j /2E 0 and E k (p) ≈ E 0 + ξm 2 k /2E 0 , respectively. The dimensionless constant ξ is smaller than one and it entirely depends on the particulars of the neutrino production process, such as the characteristics of the interaction and the nature of the other particles taking part in the process. Within the wavepacket treatment of the neutrino production process, for instance, the constant ξ can be fully expressed in terms of the dispersion width σ P of the produced neutrino's wavepacket and the corresponding group velocity of the latter [48]. By expanding these conserved energies up to the first order in m 2 /E, however, one would only end up with correction terms that are of the second order in m 2 /E in Eqs. (17) and (18) where such terms have already been discarded. Therefore, both approaches are consistent and lead to the same final result. Furthermore, in the limit Ω(φ) → 1 our results reduce to the ones obtained in previous works on neutrinos not conformally coupled [42,44].
It is also worth noting here the following. If we had not taken into account the fact that all the mass eigenstates are detected at the same time, instead of the common time integral T * given by Eq. (20), we would have inserted into each phaseΦ PW j the time integral T j obtained by using the corresponding link (19) between dt and dr. This would have led us to the same issue already discussed by various authors in the literature, which gave rise to a controversy over the presence or absence of a multiplicative factor of 2 in the oscillation length (18) [49][50][51][52][53][54][55][56]. It should be noted here that the correct and more rigorous derivation of the gravitationally induced neutrino flavor oscillation was first conducted in Ref. [49] for neutrinos not conformally coupled. The latter study has more recently been elaborated further and applied in Ref. [57] to study flavor interference patterns of neutrinos scattering off a static black hole by combining both analytic and numerical methods which allows one to include even non-radially propagating neutrinos.

III. CONFORMALLY COUPLED WAVEPACKET NEUTRINOS
The wavepacket approach to neutrino flavor oscillations in the absence of conformal coupling aims at taking into account the uncertainties associated with the production and detection processes [42]. It consists in replacing the plane waves |ν 0 j (t, x) of each mass eigenstate j by a superposition (usually chosen to be a Gaussian) of plane waves of infinitely many possible values of the momentum p j , all centered on a mean valuep j . Thus, in the wavepacket approach the normalized wavefunctions of the produced and the detected neutrinos at any spacetime point (t, x) are taken to be The Gaussian dispersion width σ p should be distinguished for produced and detected neutrinos by introducing two different widths σ pP and σ pD , respectively [42]. However, for the sake of clarity and simplicity, in what follows we only take into consideration one dispersion width σ p , to be associated with the production process. Therefore, we assume that the detected mass eigenstate state |ν k is just one of the orthonormal basis states of the Hilbert space. Note, also, that we do not consider a covariant Gaussian wavepacket, the exponential factor of which would have been exp [58]. The reason is that, doing so, would only bring corrections to our results that are of orders higher than O(m 4 j /E 2 j ), whereas we restrict ourselves here only to the first order in m 2 j /E j . As our goal here is to examine only radially propagating neutrinos, we shall consider again the static and spherically symmetric metric (10), and restrict our study to the one-dimensional version of the wavepacket (23). For radially propagating conformally coupled neutrinos within a spacetime of metric g µν , the Gaussian wavepacket (23) takes the following form in one dimension: Now, the subtlety we alluded to in the Introduction concerning the wavepacket approach resides in the two possible ways of computing the accumulated phase of the wavepacket (24) when working in curved spacetimes. It was shown in Ref. [44] that to properly implement the wavepacket approach in curved spacetimes, one needs to first locally evaluate the Gaussian integral (24) and then compute its accumulated phase at the detection point. For this reason, we shall compute the transition probability using this approach, but for the sake of completeness we shall discuss in detail in subsection III C why a second approach for dealing with the wavepackets leads one further astray when working with conformally coupled neutrinos. The first approach consists then of building the Gaussian wavepacket by first evaluating integral (24) for all values of t and r. To carry out the integration, we need to express the energyẼ j (p) in the exponential explicitly in terms of the momentump j . Since the Gaussian factor is rapidly suppressed away from the mean momentump j , we may, to a good approximation, just expandẼ j (p) in terms ofp j aroundp j to a few powers of (p j −p j ). As shown in Ref. [44], in order to obtain the term responsible for the spreading over time of the wavepacket, one needs to expandẼ j (p) at least up to the second power of (p j −p j ). However, since the resulting corrections are of the order ofm 2 j /Ẽ 3 j (p) and the resulting expressions are very involved even for the case of neutrinos not conformally coupled [44], we shall restrict ourselves here to the linear term in the expansion ofẼ j (p). The effect of the conformal coupling of neutrinos on their flavor oscillations will indeed be clearly identified and more easily confronted to experiments even at this order. Let us therefore perform the following Taylor expansion of the energy: Here,Ē j represents an average energy for the wavepacket associated to the mass eigenstate |ν j , and˜ j represents the group velocity of the wavepacket. To evaluate the parameter j (r) which is position-dependent, we first substitute the dr/dτ we extract from the second identity in Eq. (11) into Eq. (13) to find,Ẽ From this identity, we then find the following expression for Substituting this expression into Eq. (25), and then plugging the latter into the Gaussian integral (24), we arrive at the following expression for the wavepacket at any spacetime position (t, r): We need now to evaluate, using the Stodolsky prescription, the accumulated quantum mechanical phase of this wavepacket as it propagates towards the detector at r = r B . For that purpose, we need to integrate over t and r only those terms that do contribute to the quantum phase [44]. On the other hand, just like in the case of plane waves studied in Sec. II, the average energy in the phase term −iĒ j t is not the conserved quantity along the path of the wavepacket. The conserved average energy now isĒ j =Ē j /Ω(φ), as imposed by Eq. (12). Therefore, when performing the time integral to get the accumulated quantum phase we should integrate Ω(φ)Ē j by moving onlyĒ j out of the integral and leaving inside the latter the conformal factor Ω(φ). However, unlike in the plane wave case, doing so would give here two different times inside the wavepacket (28): the time T of Sec. II in the phase part of the exponential and a total time T in the decaying part of the exponential. With two different times at hand, there would be no possibility for averaging the final transition probability over time as we did in Sec. II.
There are two different possible ways of solving this issue. The first, consists of averaging the wavepacket (28) over the time t before computing the quantum mechanical phase the wavepacket would accumulate as it propagates towards the detector. The second way consists of computing first the final quantum phase of the wavepacket by integrating the phase part of expression (28) over t and r after taking into account an analog of relation (19) between dt and dr. We shall now examine both methods in detail and discuss their respective outcomes, starting with the first method.
A. Averaging over t Averaging over time by integrating the Gaussian (28) over the coordinate time t, we easily find the following timeindependent expression for the wavepacket at any radial position r: Note that this result is valid only for a time-independent conformal factor Ω(φ). To find the final value of this wavepacket at the detector, we need to compute the accumulated quantum mechanical phaseà la Stodolsky by integrating the imaginary part in the exponential while simply evaluating the real term at r = r B . To properly compute the quantum phase, we make use of Eqs. (14) and (27) to find the following expression by keeping only terms up to the first order inm 2 j /Ē j : In the second step we used Eqs. (3) and (12), and in the third step we introduced the notation K for the new integral involving the conformal factor. With this accumulated phase at the detection point, the wavepacket (29) leads to the following transition probability: Here, the algebraic oscillation length L jk osc is again given by Eq. (18), whereas j and k together with Ω(φ B ) are all evaluated at r = r B . Note that the average energy E 0 of the massless-neutrinos approximation appearing in Eq. (18) comes here again from the approximationĒ j = E 0 + ξm 2 j /2E 0 that arises from a quantum field theoretical treatment of neutrino wavepackets [48]. The dimensionless constant ξ being less than one [48], we make here and in the next subsection the approximationĒ j ≈ E 0 that restricts us to the first order in m 2 j /E 0 . Note also that we have used the fact that˜ j = j ,˜ k = k andσ 2 p = σ 2 p . The conformal invariance of j (r) can be seen by inserting Eq. (4) and the first identity in Eq. (3) into Eq. (27). The conformal invariance of σ 2 p =g µν σ µ p σ ν p stems from the fact thatg µν = Ω 2 (φ)g µν and the fact thatσ µ p = Ω(φ)σ p , as it follows from Eq. (4); leading toσ 2 p = g µνσ µ pσ ν p = σ 2 p . We clearly see now that in the limit Ω(φ) → 1, the result (31) does not reduce to what is found for the case of neutrinos not conformally coupled. This shows that this approach is not valid.

B. Averaging over T
Because of the presence of the conformal factor Ω(φ), we compute here the accumulated quantum mechanical phase of the wavepacket (28) by turning the time integral of the phase term into a spatial integral thanks to the link (19). According to our discussion in subsection II B, we should here take care again of inserting the analog of the unique time integral (20) of the most massive mass eigenstate into the quantum phase of the wavepacket of each mass eigenstate. The analog of the integral (20) for wavepackets is obtained by replacing E * (p) in Eq. (20) by the average energyĒ * of the wavepacket of the most massive state. Therefore, the expression of the wavepacket of the j-th mass eigenstate at the detection point r = r B , such that r B − r A = L, takes the following form We have used here the conformal invariance of both σ 2 p and j as discussed in the previous subsection, and we have introduced the unmeasured total coordinate time T = B A dt. Therefore, adopting again the zeroth-order approximation for the average energies,Ē j = E 0 +O(m 2 j /E 0 ), the transition probability takes, up to the first order in m 2 j /E 0 , the following form Note that, as in the case of the plane waves dealt with in Sec. II, there is no overall multiplicative factor involving Ω(φ) or its powers. Furthermore, the oscillation length L jk osc in Eq. (33) is given again by the standard expression (18). Averaging now this probability over the unmeasured time T , we easily find the following final expression for the transition probability: (34) We have introduced here the algebraic coherence length 3 : In the second step we used expression (27) for the group velocities j and k . This coherence length is nothing but the coherence length already found in Ref. [44] for neutrinos not conformally coupled. The quadratic term (L/L jk coh ) 2 in Eq. (34) is the term responsible for damping the flavor oscillations with increasing distance due to the spreading of the wavepacket.
Very important to note also, is the fact that this result shows that the conformal factor Ω(φ) does not contribute to the damping of the flavor oscillations. The conformal coupling of neutrinos affects thus only the oscillation phase. Furthermore, since the result (34) remains valid no matter how slowly the field φ varies with position, it should remain valid for a constant factor Ω(φ 0 ) achieved when the field φ 0 settles down to its minimum value φ within dense environments. Therefore, we see that even the wavepacket approach to neutrino flavor oscillations indicates that the latter are also subject to a screening mechanism. Of course, the result (34) also remains valid in the limit Ω(φ) → 1, for which we simply recover the result obtained in Ref. [44] for the case of neutrinos not conformally coupled.

C. Building the wavepacket at the detector
The second approach for dealing with wavepackets in curved spacetimes consists of building the wavepacket of each mass eigenstate from plane waves only after the latter have all reached the detector [43]. Therefore, the Gaussian integral one needs to evaluate would be the following Here, the final accumulated quantum phaseΦ PW j of each constituent plane wave is thus computed before building the wavepacket. We have now two different options again for computingΦ PW j . One may compute the quantum phaseΦ PW j either by integrating independently over t and r as done in subsection II A, or by converting the time integral into a spatial integral as done in subsection II B. In the first case, one should insert into the Gaussian integral (36) the quantum phase given by Eq. (15), whereas in the second case one should insert into the Gaussian (36) the quantum phase given by Eq. (21). We clearly see now the big issue one faces when adopting this approach. Both phases (15) and (21) would prevent one from properly performing the Gaussian integral (36). Indeed, the integration variable in that Gaussian is the momentump j in which the implicit conformal factor Ω(φ) should now be evaluated at r = r B , whereas both phases (15) and (21) contain through the energies E j (p) and E * (p) a momentum p j that is decoupled from the conformal factor Ω(φ) as the latter has already been integrated over the path of the neutrinos inside integrals T , I and J. This shows that this approach is also inadequate when applied to conformally coupled neutrinos.

IV. TWO-FLAVOR NEUTRINOS CASE
All the results we derived in the previous sections were general and apply even to the three-flavor neutrinos case. However, for the sake of clarity and simplicity, we shall apply now those results to the case of two flavors, and hence also to two mass eigenstates j, k ∈ {1, 2}. In addition, we consider here only transitions among electron neutrinos ν e and muon neutrinos ν µ , so that α, β ∈ {e, µ}. The ν e − ν τ transition probabilities are computed in exactly the same way. This case is also the most relevant in view of its application to solar neutrinos. The mixing matrix U * α j in the two-flavor case has the following form in terms of the mixing angle θ [42]: By plugging this matrix into Eq. (17), we find the following transition probability in vacuum within the plane-wave approach for a given coordinate sparation r B − r A : Here, the length J is the one given by the second integral in Eq. (15), and the algebraic oscillation length L jk osc as given by Eq. (18) has been replaced here by L 21 osc . From this expression, the survival probability can also be computed using the identity, P PW ν e →ν e = 1 − P PW ν e →ν µ . We would like to stress here that the validity of the latter identity for computing the survival probability holds only thanks to the conservation of unitarity under the conformal coupling as discussed in Secs. II and III. On the other hand, plugging the matrix (37) into Eq. (34), we find the following transition probability in vacuum within the wavepacket approach for the same coordinate distance r B − r A = L: The coherence length here is given by For a concrete numerical example, we shall plot the variation of both probabilities with the distance traveled by the neutrinos within the chameleon model for which a detailed profile inside a constant-density and compact spherical object has been derived in Ref. [59]. In what follows, we take ∆m 2 12 = 10 −4 eV 2 , E 0 = 0.5 MeV, sin 2 2θ = 1 [60] and σ x ∼ 10 −8 cm [61] for which we find L 21 coh ∼ 10 8 m. In addition, we take the scalar field's profile in vacuum outside the uniform-density star of radius R (not necessarily the solar radius) to be given by the following equation [59]: Here, φ ∞ is the value of the scalar field outside the star, m ∞ is the effective mass of the scalar field outside the star, x = r/R so that x ≥ 1, and K is a constant given by [59]: .
The value of the scalar field inside the star is denoted here by φ c and m c = βρ c /M p φ c is related to the effective mass of the scalar field inside the star [59]. These formulas are valid approximations only when a constant density is assumed for the star inside of which the scalar field is assumed to vary radially, and m ∞ R 1 [59]. On the other hand, outside the star the spacetime metric is given by the Schwarzschild solution for which A(r)B(r) = 1. Therefore, with the profile (41) for the scalar field, integral J is computed as follows In the second line, we set r A = R and we denoted by x B the ratio r B /R . In the last line, we assumed m ∞ R ≈ 0 [59], and we introduced the exponential integral function We now take the radius R to be the radius of the Sun, the average density of the star to also be of the same order as that of the Sun, ρ c ∼ 1 g/cm 3 , and the density of the environment outside the star to be that of the interplanetary medium, ρ ∞ ∼ 10 −24 g/cm 3 . Therefore, for a chameleon field with an n = 3 power-law potential, we have φ ∞ = (ρ c /ρ ∞ ) φ ∞ ∼ M p and β = 1 [62], integral J approximates to We have introduced here the coordinate distance L = r B − r A and we have kept in the second step only the leading-order terms which allowed us to discard the difference between the last two terms inside the square brackets. This shorter expression is much easier to work with. Note that for values of φ ∞ much smaller than M p , integral J would just reduce to the usual coordinate distance L, for which the transition probability would reduce to that of neutrinos not conformally coupled, and the effect of the scalar field would thus become totally screened. Therefore, the effect of the scalar field on neutrino flavor oscillations is to give rise to an effective path difference for neutrino mass eigenstates. The latter experience thus a longer path difference because of their interaction with the scalar field. Physically, this is very reminiscent of light rays propagating inside a refractive medium: two light rays experience a longer path difference when one of the rays propagates inside a medium of a higher refractive index. In fact, plugging our numerical values into the transition probability (38) for a completely screened scalar field gives whereas plugging expression (44) where L displayed in these expressions is measured in meters. The variation of the two transition probabilities P PW 0 ν e →ν µ (L) and P PW ν e →ν µ (L) with the distance L is plotted in Fig . 1 below. The effect of the conformal coupling of neutrinos to the chameleon field is thus to modify indeed the effective traveled path of the neutrinos towards the detector as illustrated by the small spatial lag acquired by the thick line representing 4 Note that considering a more compact object like a neutron star instead of the Sun has the advantage of providing us with a stronger gravitational field that would show up in a smaller radius R and a higher central density ρ c entering our formulas (41), (42) and (43). This would indeed remove the need to explore very large distances from the source to witness the effect. Fortunately, however, when considering the mean oscillation probability, as we do in Fig. 3, our physical results become much clearer graphically even when restricting ourselves to the gravitational field of the Sun. P PW ν e →ν µ (L) with respect to the thin line representing P PW 0 ν e →ν µ (L) in Fig. 1.

Transition probability
FIG. 2. The variation (with the coordinate distance L) of the transition probability P WP νe→νµ (L) given by expression (48) within the wavepacket approach is represented by the thick line. The thin line represents the variation of the transition probability P WP 0 νe→νµ (L) of neutrinos not conformally coupled as given by Eq. (47). .
Similarly, by plugging our numerical values into the transition probability (39) for a completely screened scalar field within the wavepacket approach, we find whereas plugging our numerical values as well as expression (44) for the J integral into Eq. (39), gives the following expression for the transition probability of conformally coupled neutrinos within the wavepacket approach: × cos L 10 3 + 10 6 exp − 10 7 L + 10 9 − 10 6 e 0.01 .
The variation of the transition probabilities P WP ν e →ν µ (L) and P WP 0 ν e →ν µ (L) with the distance L are plotted in Fig . 2. The exponential damping displayed in the graphs is the same as the one experienced by neutrinos not conformally coupled, and is thus not affected by the conformal coupling of the neutrinos. Only the oscillation frequency is affected. Indeed, just like in the case of neutrinos not conformally coupled, we see that the transition probability tends to 1 2 as the coordinate distance L goes to infinity. This means that at long distances flavor oscillations become washed out by the damping due to the loss of coherence between the mass eigenstates' wavepackets. The probabilities for detecting either flavors become then equal. In other words, the transition probability at long distances tends to the value expected from a statistical mixture of the flavor states.
Now, because the transition probabilities (38) and (39) are highly oscillatory and it is not possible in practice to measure those oscillations for precise values of neutrino energies E 0 and traveled distances L, we shall extract from expressions (38) and (39) the oscillation probabilities averaged over a specific distribution of L and E 0 . For that purpose, we shall follow Ref. [42] and assume, for simplicity, a normal distribution of L/E 0 , with the mean value L/E 0 and standard deviation σ L/E 0 . We assume σ L/E 0 to be proportional to L/E 0 and take for definiteness σ L/E 0 = 0.2 L/E 0 as in Ref. [42].
Starting with Eq. (38), we set J = L for the completely screened model and keep the same numerical values we chose previously, except that now we leave the energy E 0 as a variable measured in MeV. We then compute the average of the cosine function in that expression as follows: Plugging this expression of the average cosine into Eq. (38) gives the following mean oscillation probability for the completely screened scalar field: On the other hand, to extract the mean oscillation probability for neutrinos coupled to the chameleon field we first expand J as given by Eq. (44) up to the first order in L and get J ≈ eL. Therefore, the average oscillation probability in this case reads P PW ν e →ν µ = (51) Both transition probabilities as functions of the mean ratio L/E 0 are plotted in Fig. 3 below. We much more clearly see from this plot the effect of the scalar field on the flavor oscillations. The spatial lag between the two oscillation probabilities can distinctly be read off from the graph.

V. THE CASE OF NON-UNIVERSAL CONFORMAL COUPLINGS
In this section we describe how our previous results get modified when the conformal coupling of neutrinos is not universal. We say that the conformal coupling is non-universal when the scalar field φ and/or the conformal factor Ω(φ) is different for different mass eigenstates of the neutrinos. We denote here the two possibilities at once by simply adding a subscript j; denoting thus by Ω j (φ) the conformal factor that interacts with the j-th mass eigenstate.

A. Plane-wave neutrinos
When we dealt with neutrino plane waves in Sec. II, we had the option to either integrate independently over time and space, or to turn the time integral into a spatial integral by taking care of assigning to all the mass eigenstates the same arrival time T * . We argued that although the first option is not as rigorous as the second, it does nevertheless lead to the same final transition probability as the second. However, when allowing for a non-universal coupling we immediately see that the first option becomes totally inapplicable. In fact, the first option would assign to any j-th mass eigenstate a quantum phase that has the form It is clear that when substituting this phase into the expression (9) for the transition probability, no single time T would be available to be integrated over to extract an average transition probability.
If, on the other hand, we use the second option then it becomes possible to compute the transition probability by using again the common time given by Eq. (20). However, such a time, corresponding to the most massive eigenstate, takes now the following form: Replacing the time integral T j by this time T * in the calculation (53) of phases, the resulting quantum mechanical phase a j-th mass eigenstate acquires is given by the following expression: Therefore, inserting into expression (9) of the transition probability the phase differenceΦ PW jk one obtains by subtracting from Eq. (54) the corresponding expression for the k-th mass eigenstate, and then using E j (p) ≈ E 0 + ξm 2 j /2E 0 and E k (p) ≈ E 0 + ξm 2 k /2E 0 , the transition probability takes the following form We introduced here the notation ∆I jk ≡ I j − I k . The oscillation length L jk osc is still as given by Eq. (18). Note that we have assumed here that the mixing matrix U * α j is the same for conformally coupled neutrinos even in the non-universal conformal coupling case (see Appendix A for the case where this assumption is relaxed). In contrast to the universal coupling scenario, however, we see that the parameter ξ appeared here explicitly in the transition probability even though we limited ourselves to the first order in m 2 j /E 0 . From the result (55) we find the following probability for an electron neutrino to be detected as a muon neutrino in the simple case of two-flavor neutrinos: This result neatly shows how a non-universal conformal coupling would affect the transition probability. It manifests itself in the difference between the ways each mass eigenstate experiences the spacetime path along which it travels. Because of the very last term inside the parentheses, we see that the correction becomes more significant for high-energy neutrinos. When spacetime paths are identically experienced by all the mass eigenstates, i.e., when J 2 = J 1 and I 2 = I 1 = I * , we recover the result (38) of the universal coupling scenario.

B. Wavepacket neutrinos
When dealing with neutrino wavepackets in the nonuniversal conformal coupling, the calculation proceeds along the same steps that led us to Eq. (32) and then to Eq. (33). We compute here the accumulated quantum mechanical phase of the wavepacket (28) by using the same recipe we used for computing the phase (54) of each mass eigenstate. The result has then the same form as in Eq. (54), except that E j (p) is replaced by the average energyĒ j and E * (p) is replaced by the average energyĒ * of the wavepacket of the most massive state. Consequently, the expression of the wavepacket of the j-th mass eigenstate at the detection point r = r B , such that r B − r A = L, takes the following form Therefore, adopting again the first-order approximation for the average energiesĒ j = E 0 + ξm 2 j /E 0 , the transition probability takes, up to the first order in m 2 j /E 0 , the following form Here, the oscillation length L jk osc is also given by the standard expression (18). Averaging now this probability over the unmeasured coordinate time T , we easily find the following final expression for the transition probability: The algebraic coherence length L jk coh is as given by Eq. (40). Note that the non-universality of the conformal coupling does not affect the coherence length either. From this result, we compute the transition probability for the two-flavor neutrinos case to be, This result reduces, of course, to expression (39) when the conformal coupling becomes universal.

VI. CONCLUSION AND DISCUSSION
We have examined the effect a conformal coupling of neutrinos would have on their flavor oscillations when they are propagating in vacuum. We adopted a general static and spherically symmetric metric for the curved spacetime and derived the transition probabilities in vacuum both within the plane-wave approach and within the wavepacket approach. We showed that in this case there are two different ways of dealing with time in the plane-wave approach: Either by averaging the transition probability over the unmeasured time as usually done in the literature for neutrinos not conformally coupled, or by converting the time integral into a spatial integral and taking into account the difference of the times of arrival of the different mass eigenstates. We showed that both methods lead to the same correct oscillation length. However, we showed that the first method is not applicable when dealing with wavepackets of conformally coupled neutrinos.
We have discussed in great detail the different possible ways of implementing the wavepacket approach and highlighted the additional subtleties the latter conceals due to the presence of the position-dependent scalar field φ along the path of the neutrinos. For simplicity, we studied only the case of a time-independent scalar field φ since the calculations would become unwieldy when the field φ depends on time. Then we examined in detail the possibility of implementing the approach that consists of evaluating the wavepacket at the detector and we showed that the issues it gives rise to are more serious than the ones that already arise in the absence of any conformal coupling. We numerically examined the simple case of two-flavor neutrinos and plotted the variation with distance of the transition probabilities within the chameleon model using both the plane-wave and wavepacket approaches.
Although screening models make it possible to violate the equivalence principle even when their scalar field couples universally to matter, we have nevertheless explored in this work the additional possibility of having a non-universal conformal coupling in those models. A non-universal coupling offers, indeed, an additional way of violating the equivalence principle. We derived the general formulas for the probability transitions for such a scenario. Our formulas showed that the violation of the equivalence principle through these flavor oscillations would be more important for high-energy neutrinos.
Finally, we would like to stress here the important fact that all the transition probabilities derived here never exceed unity, both with a universal and a non-universal coupling scenario. This preservation of unitarity is in fact what gives meaning to the flavor survival probability and allows one to compute the latter. We pointed out that this important conclusion is due to the fact that a conformal transformation on Dirac spinors preserves unitarity.