Damped Neutrino Oscillations in a Conformal Coupling Model

Flavor transitions of Neutrinos with a nonstandard interaction are studied. We consider a model with a scalar field conformally coupled to matter and neutrinos. This interaction alters the neutrino effective mass and also its wave-function leading to a damping factor, causing deficits in the probability densities and affecting the oscillation phase. As the matter density determines the behavior of the scalar field, we also have an indirect matter density effect on the flavor conversion. We explain our results in the context of screening models and study the neutrinos oscillation and the deficit in the total flux of electronneutrinos produced in the Sun through the decay process, and confront our results with observational data.


Introduction
One of the interesting subjects in particle physics and also in cosmology is the physics of neutrinos which has received a lot of attention especially in the beyond standard model theories. In the standard model of particles, neutrinos were assumed to be massless. Massive neutrinos were first proposed theoretically in [1,2] and after some decades finds severe evidence by detection of deficits in the number of electron neutrinos received from the Sun [3,4], and in the number of atmospheric muon neutrinos [5,6]. The discrepancy between the theoretically expected number of neutrinos and the observation may be explained by considering the flavor states as mixing of mass state neutrinos. This gives rise to flavor change via oscillation during their travel (in the atmosphere), and also adiabatic flavor conversion (as in the Sun) through the MSW (Mikheyev-Smirnov-Wolfenstein) effect describing the influence of local electron density on electron neutrino density via weak interaction [7]. Apart from the standard weak interaction, one may consider nonstandard neutrino interactions (not included in the standard model) [8] and investigate their influences on neutrino oscillations [9]. Non-standard interaction of neutrinos and exotic fields has also attracted many attentions in the cosmology [10][11][12][13][14]. These additional fields may be the dark energy describing the present cosmic acceleration of the Universe. Inspired by the similarity of the neutrino mass squared difference the scale of dark energy, mass varying neutrinos (whose mass depends on a light scalar field) were introduced in [15][16][17] and used in [18,19] to discuss the solar neutrino flavor survival probability. The dependence of the neutrino mass on the environment due to the interactions with dark sectors was also studied in [20], showing that the interaction, by affecting the neutrino mass, alters drastically its oscillation.
In [21][22][23], a model for the Universe cosmic acceleration [24] based on quintessence-neutrino interaction through a conformal coupling was proposed such that when the massive neutrinos became non-relativistic, by activating the quintessence, ignited the Universe acceleration. In this class of scalar-tensor dark energy models, the quintessence is interacting with matter (including neutrinos)through a conformal coupling [25][26][27]. This coupling may give rise to the screening effect as was studied in the chameleon [28][29][30][31]and the symmetron [32,33] models. In these screening models, the behavior of the scalar field is specified by the matter density such that in a dense region they are screened. In this paper, we aim to study the effect of such a conformal coupling on the densities and oscillations of neutrinos. As the behavior of the scalar field is specified by the matter density, we expect to encounter an MSW-like effect but caused by a non-standard interaction i.e. the conformal coupling. The interaction with the scalar field by modifying the neutrino effective mass and its density provides a mixed situation [34] in which the neutrino deficit may be related to both the neutrino oscillation and neutrino decay. This paper is organized as follows: In sec.2 we discuss the ingredients of our model starting with a conformal transformation applied on the metric which rescales relevant Dirac action governing the motion of neutrinos. Corrections to the neutrino mass and wave-function related to the curved spacetime in the presence of a scalar field are explicitly calculated in this section. We then give an analytical solution to the mass-varying Dirac equation in subsec.2.1 and show how the scalar-matter coupling function will change the neutrino quantum mechanical phase. The damped two-flavor transition is the subject of the subsection2.2 giving the probabilities for both chameleon and symmetron cases. After that, we will extend damped neutrino transitions to three flavors as well as discussing the possibility of neutrino decays and deficit in the total probability in 2.3. Interactions with matter such as the Sun, for instance, have a dramatic effect on the flavor conversion, MSW (Mikheyev-Smirnov-Wolfenstein) effect, presented in the subsection2.4. We will see that the LMA (large mixing angle) solutions have the most plausible effect on the neutrino decay. Finally, we summarize and present our results in sec.3 through some numerical examples. A comparison of our results to the SNO (Sudbury Neutrino Observatory) [35] and also to the Borexino [36,37] data for the MSW-LMA survival probability will be done.

Conformal coupling and neutrino oscillation
Our model is specified by the following action including a scalar field conformally coupled to matter where g is the determinant of the metric, g µν , M p is the reduced Planck mass, R is the Ricci scalar and L m is the Lagrangian density of other components generally denoted by Ψ i .g µν is related to the metric g µν , by [38][39][40] The conformal factor A(φ) is a function of φ. Various screening models are resulted by different choices of the coupling and potential functions, e.g. the chameleon model corresponds generally to the choice is a field-dependent coupling parameter. A simple case corresponds to a constant value of β ∼ 1. In the symmetron model a quadratic coupling function such as A(φ) ≡ 1 + φ 2 (r) 2M 2 is chosen to respect Z 2 symmetry. In the Appendix A, we will review these two screening mechanisms and point out the required relations for our discussion. If we want to respect the weak equivalence principle (WEP), we must take a universal A(φ) [28], otherwise, we may have different A i (φ) corresponding to different ψ i .
To compute the probability of neutrino flavor transition, we study the neutrino equation of motion. The neutrino action is To obtain the covariant derivative we proceed as follows: The vierbein or tetrad, µ a , arẽ The Dirac γ-matrices and the Christoffel connections Γ λ µν are given bỹ respectively. The spin connection is in which ω ab µ = − bν (∂ µ a ν − Γ λ µν a λ ). So the covariant derivative is The commutation of γ-matrices, is unchanged under the conformal transformation: The second term in (8) can be written in a more familiar form bν a Therefore the covariant derivativeD µ can finally be written as leading toγ One can use this relation to show that provided that and In the following we fix g µν as the Minkowski metric η µν and the covariant derivative D µ simply becomes the partial derivative ∂ µ .

Neutrino oscillation
Neutrino oscillation or mixing can be explained in terms of the relationship between flavor eigenstates (ν e , ν µ , ν τ ) and mass eigenstates (ν 1 , ν 2 , ν 3 ). This implies that we can write each flavor eigenstate as a superposition of three mass eigenstates. Mixing is given by where  are mass and flavor eigenstates, respectively, and U which is a unitary mixing matrix, called PMNS (Pontecorvo-Maki-Nakagawa-Sakata) matrix, can be written as Dirac Lagrangian density describing both right-handed and left-handed neutrinos including the mass and dynamical terms is given by [41] where m αβ 's are components of a 3 × 3 mass matrix in the flavor basis, ν α(R,L) are 2-spinors implying right-handed and left-handed neutrinos of flavor α. In addition, σ µ R ≡ (σ 0 , σ 1 , σ 2 , σ 3 ) and σ µ L ≡ (σ 0 , −σ 1 , −σ 2 , −σ 3 ) are Pauli matrices for right-handed and left-handed spinors, respectively. We also note that all primes in the formulas refer to the rescaled relations. Now, by varying this Lagrangian with respect to the spinors, equations of motion can be obtained. For the first one, we have and the second equation is given by Generally, the left-handed and right-handed spinors can be defined by Substituting these relations into Eqs. (19) and (20), we obtain Notice that derivatives in relations above are all taken with respect to the radial coordinate r, because we have only considered the radial propagation in our model to solve Eqs. (23) and (24). For neutrinos masses much less than their energies, the second term in the square brackets in Eq.(24)is negligible with respect to the first one. By applying this approximation, we obtain and by replacing this into Eq.(23), we have where m βγ = U * βi m i and m * αγ = U αi m i . This leads us to the following equation By using f i (r) = U αi f α (r) and also Eq. (27) we obtain This equation can be easily solved by Therefore, according to Eq.(21), the normalized space-time part of the neutrino wavefuction denoted by Ψ(r, t) satisfying the initial condition is given by (see also [42][43][44][45][46] where As we will see depending on the model, D may be a damping or an enhancing factor. We will note it by D factor. is the phase of the oscillation. Note that in (31) and (32), we assume generally different A for different neutrinos. By A 0i we mean the value of A i at the initial point (r 0 , t 0 ). Notice that the initial value of the function F i (r, t) is clearly equal to one, i.e. F i (r 0 , t 0 ) = 1. In what follows, the above results will be used to calculate the oscillation probability.
where we have used F i (r 0 , t 0 ) = 1. Choosing the initial condition as Ψ 1 (r 0 , t 0 ) = Ψ 2 (r 0 , t 0 ) =: Ψ 0 , we find where we have used |ν e = cos θ |ν 1 + sin θ |ν 2 . At the moment, we pay attention to the probability of neutrino oscillation for the case of two-flavor neutrinos, e.g. ν e and ν µ . As we know, the relation between mass and flavor eigenstates can be described by a 2 × 2 mixing matrix as follows Using (36), the neutrino state (33) can be written in terms of flavor eigenstates as The transition probability ratio is given by and by substituting from Eq.(37) we have where Φ 12 (r) = ϕ 1 − ϕ 2 is the phase difference between different neutrino mass eigenstates. The survival probability ratio, on the other hand, can be obtained by multiplying (37) by ν e | from left-hand-side, so we have and doing some manipulation leads to These results are similar to the results obtained by Refs. [47,48]. The conformal coupling beside the influence on the oscillation phase through Φ 12 has a damping or enhancing effect via D ij : which is not generally equal to one. For D ij = 1, we have only oscillatory wavefunctions and the sum (42) equals one. Besides the neutrino flavor oscillation we may have neutrinos production or annihilation via the scalar field interaction. Note that if the conformal coupling is the same for all the neutrinos, we have D ij = Hence where . If both mass eigenstates are coupled with a same coupling to the scalar field φ, i.e. β 1 = β 2 = β, the probabilities (44) will reduce to To obtain the values of the above probabilities, we can use the solution to the chameleon equation of motion, which is solved numerically in subsection A.1. According to relations (45), not only the scalar field changes the oscillation phase but multiplies the probabilities by a scalar field dependent coefficient (for β 1 = β 2 this coefficient becomes relevant for the relative number of flavors to each other). The sum of the probabilities is not one unless φ(r) = φ(r 0 ). Hence generally implying neutrino-scalar interaction, affecting the neutrino density.
Similarly for the symmetron model, the factor D is Using these relations leads us to the following probabilities We can also obtain the sum of probabilities above which is independent of mixing angle θ and the phase shift Φ 12 , Provided that the field φ(r) > φ(r 0 ), the sum of the probabilities is smaller than unity and D may be interpreted as a damping factor and conversely for φ(r) < φ(r 0 ) the neutrino density is enhanced.

Three-flavor oscillations
We extend our results by taking three flavors into the account. As before, we first have to write the most general state of the neutrinos. The evolved neutrino state corresponding to the flavor α can be suggested as follows where F i (r, t) is a function which consists of D and phase factors. Note that F i (r 0 , t 0 ) = 1, which explicitly implies that the initial condition can be satisfied by (50). In order to obtain the different probabilities for this case, we shall replace the mass eigenstates by flavor eigenstates. As the two-flavor case, the mixing between various flavors of neutrinos can be described by a unitary mixing matrix U as follows or equivalently Using relations (51) and (52), the state (50) can be written as Assuming that at (r 0 , t 0 ) we have a specific flavor neutrino e.g. ν α , requires that Ψ i (r 0 , t 0 ) = Ψ 0 for all i's. The probability amplitude ratio (α → β) is then Therefore, the probability formula is given by Now, we can calculate this formula in more details by using the explicit forms of the functions F i (r, t) Eq.(30), thus we have where and by doing some manipulation, the relation (56) can be rewritten as where the two first non-oscillatory terms belong to the damped(enhancing)neutrino mixing. From another standpoint, the last term includes imaginary sector of the mixing matrix, so this term might correspond to the CP-violation.
To calculate probabilities P ee , P eµ , etc., we have to apply the general form of the mixing matrix [49] where c ij ≡ cos θ ij , s ij ≡ sin θ ij and δ is the CP-violating phase. As an example, ν e -survival probability is given by where Σ ij (r) is defined below Eq. (44). As can be seen, the probabilities for the three-flavor case are so tedious to calculate directly for the matrix elements of (58) and they are beyond our scope in this paper. Instead, we are going to investigate a simpler case in the next subsection.

Neutrino decay: stable ν 1 and unstable ν 2
We assume that the lightest neutrino mass eigenstate, ν 1 , is stable during neutrino propagation, i.e. D 11 = 1, whereas decay of ν 2 as well as mixing among three neutrino eigenstates have effect on the deficit in the solar neutrino flux [50]. From standpoint of our model, this statement means that only one of the mass eigenstates is coupled and so affected by the scalar field. For this case, we suppose that the U e3 element of the matrix (58) is approximately negligible. This is because based on the CHOOZ experiment results one can take two possibilities for this element: U 2 e3 < 0.03 or U 2 e3 > 0.97 [51]. The solar neutrino data [52], however, excludes the large values of U e3 , so we take the first and much smaller inequality of U 2 e3 . Therefore, the PMNS mixing matrix (58) reduces to Assuming |ν e as the initial state and using the unitarity of the above matrix, we calculate different probabilities from Eq.(57) as follows where D ≡ D 22 is the D factor corresponding to ν 2 -decay. The observed neutrino rates in the SNO experiment from the CC (charged current), NC (neutral current) and ES (elastic scattering) interactions can be combined to place constraints on the separate ϕ(ν e ) and ϕ(ν µ ) + ϕ(ν τ ) fluxes, which are respectively related to P ee and P eµ + P eτ . Both of these probabilities and also sum of them, i.e. P tot = c 2 12 + s 2 12 D, are independent from atmospheric mixing angle θ 23 . On the other hand, the deficit in the rate of the solar neutrinos in this model on the Earth, i.e. δP αβ ≡ P where the above subscript X may refer to any invisible (or unknown) particle(s), e.g. a scalar field. As it can be seen, δP eX depends on the mixing angle θ 12 as well as the D factor. For increasing functions φ i (r) we have D ≤ 1 and so there is a possibility for electron-neutrinos travelling from the source to the detector to decay to φ particles. For this case, we also obtain the discrepancy between undamped and damped probabilities for three cases separately. Using relations above, we have For example if we are interested in oscillation of neutrinos produced inside the Sun, for the case δP ee > 0, we would have cos(Φ 12 ) > − tan 2 θ 12 , where we have assumed that √ D 1 outside the Sun, since both chameleon and symmetron fields are nearly constant deep in space. To have positive δP eµ and δP eτ , however, the trivial condition cos(Φ 12 ) < 1 has to be satisfied. Therefore, the resulting range for the phase difference is − tan 2 θ 12 < cos(Φ 12 ) < 1.

Flavor conversion in matter: The MSW effect
The weak interactions of neutrinos in matter modify the flavor conversion relative to the cases of propagation in the vacuum, as predicted decades ago by Mikheyev, Smirnov and Wolfenstein dubbed the MSW effect [53]. After solar neutrino production in the core of the Sun and during their travel inside the Sun, they scatter forwardly from electrons until they leave the Sun and propagate through the vacuum to the Earth and finally detected e.g. in the SNO detectors. Therefore, we should take this effect into account. We start with the following rescaled Hamiltonian of the system where A is a dimensionless parameter originated from ν e − e − scattering in matter which is defined by where G F is the weak Fermi constant and n e is the number density of electrons in the bulk of the Sun. Assuming A = cos 2θ, leads us to the resonance point indicating maximal mixing in matter. The number density at the resonance point is then given by According to Eq.(64) in the dilute regions, parameter A takes a small value so that according to the Hamiltonian, Eq.(63), this implies an ordinary mixing. In the dense regions, on the other hand, A is large enough to dominate in the diagonal elements of the Hamiltonian denoting mixing suppression in such regions, e.g. in the solar core, where neutrinos are born [42].
Remember that the rescaled mass-squared splitting ∆m 2 (r) depends explicitly on the scalar field values, see Eq. (14). The chameleon screening mechanism is specified by its varying mass so that in the dense regions, where parameter n e is considerable, this scalar field has a large mass hiding it from detectors. From properties of the chameleon (see the appendix), we see that chameleon takes small values φ min,in M p in such regions and, consequently, ∆m 2 −→ ∆m 2 . In the dilute regions, however, where the chameleon has a small mass and the parameter A becomes small, the chameleon acquires larger values φ min,out ∼ M p which results in ∆m 2 = exp(2βφ/M p )∆m 2 .
Another screening mechanism, symmetron model, is based on spontaneously symmetry breaking discussed in Appendix A.2. The Z 2 -symmetry will be broken in the low density environments which yields φ min,out = 0, see Eq.(A.15), and parameter A is small, so neutrinos experience ordinary oscillations. The restoration of the Z 2 -symmetry, on the other side, imposes φ min,in −→ 0. In such a region, A becomes dominant and, consequently, gives mixing suppression [42].
The MSW flavor conversion in the Sun can be considered as a level crossing (or resonance) at which the most flavor change occurs when neutrinos cross this point. The standard expression for "jumping"probability between ν 1 and ν 2 inside the Sun and at the resonance point is given by where the parameter α in the exponent is equal to π∆m 2 Using Eq.(66), we predict that the LMA (large mixing angle) solutions of neutrino oscillations can have the most influence on the decay process. As we know, ν e is produced mostly as the mass eigenstate ν 2 in the solar core [54]. Figure 1 shows the behavior of P J as a function of the mixing angle θ so that for the large values of this angle P J is too small, then ν 2 does not jump to ν 1 at the level-crossing. This means that the possibility of ν 2 -decay survives. On the other side, P J is appreciable for small mixing angles and, hence, ν 2 jumps to ν 1 which has been assumed to be the stable mass eigenstate. This stability destroys the neutrino decay possibility. The LMA solution is therefore has the strongest influence on the neutrino decay.
To discuss MSW effect on the probabilities [50,54] in the Sun which are a mixture of flavor oscillation and neutrino decay, we proceed as follows. As discussed in subsection 2.3.1, ν 2 -instability with mixing together cause the solar neutrino problem, whereas the mass eigenstate ν 1 remains stable, thus, P 1 and P 2 .D are respectively the probabilities of detecting ν 1 and ν 2 on the Earth. The probability P 1 (P 2 ), defined as the probability of ν e −→ ν 1 (ν 2 ), can be written as where θ m is the mixing angle in the matter defined by tan 2θ m = sin 2θ cos 2θ−A [55]. Then using the unitary mixing matrix in the subsection 2.3.1, the probability formulas are given by P ee = c 2 12 P 1 + D s 2 12 P 2 P eµ = c 2 23 s 2 12 P 1 + D c 2 23 c 2 12 P 2 P eτ = s 2 23 s 2 12 P 1 + D s 2 23 c 2 12 P 2 , where D is the D factor inside the Sun and we have also considered D vac 1. To see the damping (decay) behavior only, the phase parts in the probabilities have been ignored.

Results, discussions and conclusion
We have studied a scenario about damped neutrino oscillations in a curved space-time which consists of a scalar field conformally coupled to other ingredients such as matter and neutrinos (see Eq. (2)).
To derive the oscillation probabilities, we studied the behavior of the Dirac equation under the conformal transformation which reduces the model into the flat space-time with a re-scaled wavefunction and a coordinate dependent mass (see Eqs. (14) and (15)). By solving the mass-varying Dirac equation in the flat spacetime we derived mass-eigenstates of the neutrinos (see Eq. (30)). The presence of factor D ij (see (31)) and the deficiency in the total probability can be interpreted as the interaction of the neutrino with the scalar field allowing them to convert to each other. To be more specific we considered solar neutrinos and consider two examples: the chameleon and the symmetron models, which we review briefly and point out the required relations for our discussion in the appendix. In these models the effective masses depend explicitly on the local matter density, see Eqs.(A.7) and (A.16), leading to the screening effect in a dense area.
The dynamics of the scalar field is determined with a good approximation by the matter density. We first expand the scalar field around its background value, and obtain an equation for the fluctuation(see (A.9)). Inside the Sun, the density ( ρ (R)) is shown in 8, and the scalar field equation is numerically depicted in Figs.9 and 10 for various values of the coupling parameter β for the chameleon, and is depicted in 11 for the symmetron.
As can be seen from (32), the chameleon scalar field affects the oscillation phase. We then found the formulas for damped transition and survival probabilities and also obtained a value for violation in the total probability conservation (see Eq.(61)) for both two-flavor and three-flavor neutrino oscillations. Finally, we studied the non-oscillatory effects of neutrino forward elastic scattering from electrons in matter, MSW effect, in subsection 2.4. As we concluded, the LMA solutions of neutrino oscillations are the best solutions describing the decay process in this model. Now let us illustrate our results via some numerical examples. In Fig.2, we plot the survival probability P ee as a function of neutrino energy for the two-flavor case to see the effects of the conformal coupling on the oscillation phase. It should be noticed that we have picked the numerical values n = 1, M 2.08keV [29], the mass-squared splitting ∆m 2 = 7.4 × 10 −5 eV 2 and the mixing angle tan 2 θ = 0.41 for the LMA solution of neutrino oscillations [56]. We have also picked three values for coupling strength β in the chameleon model, panel 2a. In the chameleon model proposed by Khoury and Weltman [57], scalar field is coupled to matter with gravitational-strength, i.e. β ∼ O(1) (red curve) in agreement with expectations of the string theory [58], or smaller, i.e. β 1. Setting β ∼ 0 reproduces the no-chameleon figure (black curve). On the other hand, as the satellite experiments [59][60][61] have proposed, they are unable to put an upper bound on β [58], so we plotted a figure for β ∼ O(10 2 ) (cyan curve). The resulting profile for the survival probability is sensitive to β such that this probability will be suppressed with intense oscillations when β increases. Also, the probability P ee and its oscillating behavior affected by the symmetron mechanism is presented in the panel 2b for three values of mass-scale parameter M 10 −4 M p . This constraint on M is imposed by local tests of gravity [32,62]. Note that µ ∼ 10 −12 eV [42], λ ∼ 10 −50 (λ 10 −96 [32]).    [29]. P ee is depicted for various orders of magnitude of coupling parameter to compare the probability amplitude and phase of oscillations for different β's. Green curve is plotted ignoring the effects of the chameleon. As can be seen, the probability amplitude is suppressed with rapid oscillations for β 10 2 (cyan curve) such that it has intersection with none of experimental values for survival probability measured by SNO (blue arrow) and Borexino (red arrow), for instance. These values measured by SNO (ϕ tot = 4.94 +0. 21 −0.21 (stat) +0.38 −0.34 (syst) and ϕ(ν e ) = 1.68 +0.06 −0.06 (stat) +0.08 −0.09 (syst)) [35] and Borexino [36] correspond numerically to P SNO ee = ϕ(ν e )/ϕ tot 0.34 and P Borexino ee = 0.29 ± 0.10 within 1σ.   The behavior of the electron-neutrino survival probability is also shown as a function of chameleon-matter coupling β in Fig.5. Note that ∆m 2 = 7.4 × 10 −5 eV 2 and tan 2 θ = 0.41. The envelope function (green) shows how Pee gradually drops when β increases, due to the D factor. Fig.5 shows a damped oscillation which its damped part illustrated by the upper envelope curve (green curve) is trivial because of the D factor behind the oscillatory term which is responsible for its oscillation. As can be easily seen, P ee has rapid oscillations for increasing β which is a clear sign of the effect of the chameleon on the phase. Note that we have used the LMA mass-squared splitting and mixing angle, as mentioned before, and this figure has been plotted for solar and high energy neutrinos as shown in the legends.
Solar neutrino problem can be possibly solved by considering decay processes as well as mixing. There are generally two kinds of decays depending on whether the final state particles are only invisible, such as φ-particles, sterile neutrinos (generally non-active neutrino flavors), Majorons, etc., or they include visible particles too, e.g. active neutrino flavors [63]. An applied fit to BS05(OP) data of solar neutrinos [64] can lead us to the zeroth-order measurements 0.292 +0.067 −0.039 and 0.12 +0.14 −0.23 for the electron-neutrino survival probability and the conversion probability to unknown states respectively [65]. The latter value of the conversion probability (neutrino decay to unknown states) can be interpreted as φ-particles.
To clarify our model, we take a numerical example for δP eX . We need first to look for the numerical values of the field which minimize the effective potential for background matter density on the cosmological scales, i.e. ρ 0 ∼ 10 −24 g.cm −3 [66]. This means that we first set V eff,φ = 0 (see Eq.(A.6)), then we have . Setting parameters n = 1, β 1 and M = 2.08 keV [29] for the chameleon field gives us φ 0 ≈ 0.609074M p . As mentioned earlier, our choice of coupling parameter β is consistent with gravitation-strength interactions [67]. Using all these and also the numerical values for the symmetron field, i.e. µ = 10 −12 , λ = 10 −50 and M 10 −4 M p , we plot Figure 6 which shows the amount of the discrepancy in the damped total probability from unity for two cases of the chameleon and symmetron scalar fields. In both panels, the amount of the discrepancy grows to a constant value at R = 1 and remains constant outside the body to the Earth. The left panel of Fig.6 shows that the numerical behavior outside the sun agrees well with the numerical value presented above. The discrepancy in total probability approaches to the value δP eX 0.118. This value is much larger than that for the symmetron field depicted in the right panel of Fig.6. This difference may refer to their different coupling functions.  We have also depicted the P ee -experimental values from the Borexino data [37] of pp, 7 Be, pep and 8 B fluxes. Gray band in both panels is the best theoretical prediction of P ee (within ±1σ) according to MSW-LMA solution [37]. We guess that the best fit for this curve can be written as To draw these figures, we have assumed that ∆m 2 2.49 × 10 −4 eV 2 , tan 2 θ m 0.43 and √ 2G F n e 2.27 × 10 −7 eV 2 M eV [55]. Black, red, and cyan curves have been plotted according to the present model of MSW effect which are approximately consistent with Borexino MSW-LMA prediction (gray band). Finally, as can be seen, our model is in agreement with Borexino experimental data in low energy range (pp, 7 Be and pep) and also in high energy range ( 8 B), where changing flavor for the latter range is caused mostly by matter effects in the bulk of the Sun [37].

A Appendix: Screening mechanism
In this section, we provide a brief review of two screening models: the chameleon and the symmetron models and study their equations of motion in a nearly flat static spherically symmetric space-time.

A.1 Chameleon mechanism
This model is specified by a runaway power-law (continuously decreasing) potential of the form where n is a positive number and M is parameter of mass scale. The main feature of the chameleon scalar field is that its effective potential depends explicitly on the matter density. By varying the action (1) with respect to the field, we obtain the following equation of motion whereT µν ≡ (−2/ √g )δL m /δg µν is the energy-momentum tensor which is conserved in the Jordan framẽ From equation of statep = ωρ we know the relationship between matter density in Einstein and Jordan frames [68] where the right-hand-side of this equation can be written as the derivative of the effective potential V eff (φ) = V (φ) + A(φ)ρ with respect to the field φ where ρ is the matter density. The minimum of the potential is determined by using the equation V eff,φ (φ min ) = 0 which leads to Because the cosmological and local gravity experiments impose the condition βφ Mp 1, we assume that e βφ Mp ≈ 1. The effective mass of the field is given by where we have again used the assumption e where the first condition is for non-singularity of the scalar field at the center of spherically symmetric body, while the second implies that the field converges to a constant at infinity. The solution to the Eq.(A.5) can be obtained by expanding the field about its background as φ(r) = φ 0 + δφ up to linear order, where φ 0 is the uniform background value and δφ is the perturbation induced by a spherically symmetric body like the Sun whose radius is R . Therefore the field equation turns into where the solution to this equation in dilute regions outside the body (R := r R > 1) with constant density ρ 0 ρ is given by where δφ in (1) is the field value at the surface, i.e. R = 1, coming from the continuity condition andρ is the average density of the body. Adding this solution by the background value, we obtain the following solution for outside the body As can be seen from this, the scalar field induced by a celestial object, e.g. the Sun, is too small in large distances such that we can ignore its effects on another object in the solar system scale.
For inside the object, however, the density distribution of a spherical body like the Sun is a function of its fractional radius R ≡ r/R , as depicted in Fig.8.  Figure 8: The density distribution function of the Sun in terms of the fractional radius R ≡ r R . This figure has been drawn by BS2005 data [69].
With this density function, by solving numerically the equation (A.9) we obtain an interpolating-function with a gravitation-strength coupling β ∼ 1, as shown in Fig.9. In this figure, both the resulting perturbation and the whole field inside the Sun are depicted.  The chameleon field profile for inside the Sun (R ≤ 1) with a R-dependent mass density distribution function. Note that we have assumed that β ∼ 1, n = 1, ρ 0 10 −11 eV 4 and M 2keV [29].
We note that the resulting field profile can be sensitive to the change of the coupling β, see Fig.10. This figure shows the behavior of the chameleon scalar field. We have used the value of δφ in (R) at the surface of the Sun. The effects of the coupling parameter β on the chameleon field for four different values of β ∈ {0.1, 1, 10, 100} are shown implying that the chameleon tends to smaller asymptotic values when β grows. We also note that the allowed range of the field becomes smaller when β reduces, hence, the field is being constant for β 1.

A.2 Symmetron mechanism
In the symmetron screening mechanism, the screening is realized by symmetry restoration in sufficiently high-density regions in which φ = 0 where the symmetron-matter coupling tends to zero. In low-density regions, the Z 2 -symmetry is spontaneously broken and < φ > = 0. An example of a Z 2 symmetric coupling function and potential is: where M and µ are two parameters of mass scale and λ is a dimensionless parameter. The equation of motion of the scalar field in a static spherically symmetric background is given by where the effective potential up to a constant is as follows where ρ c ≡ µ 2 M 2 is the critical density. The breaking or restoration of Z 2 -symmetry depends on whether the matter density is smaller or larger than the critical density. In the dilute regions the symmetry is spontaneously broken and the scalar field acquires a VEV (vacuum expectation value) For the Sun, the scalar field effective mass m 2 min ≡ d 2 V eff dφ 2 | φ min is then 16) where as before ρ 0 is the background matter density andρ is the solar average density and R ≡ r R . By expanding the scalar field around φ min and keeping only the leading term we obtain the equation of motion inside and outside a body like the Sun as follows where we have used the approximation V eff (φ) ρ M 2 φ and φ 0 = φ min,out is the field value at the edge of the Sun [70]. The solution to the equation is A is a constant to be determined by continuity conditions at the boundary. For the outside, the equation of motion is given by where B is a constant. To specify the constants A and B, we should use the continuity of δφ(R) and its first derivative at the R = 1, which leads us to .
(A. 21) By assuming m in R 1 and m out R 1 and after some manipulation , the solution is derived as In Fig. 11, we have plotted the symmetron scalar field as a function of the fractional radius R. Note that we have picked the numerical values M 10 −4 M p , µ = 10 −18 eV, λ ∼ 10 −50 for model parameters. As the chameleon case, the symmetron tends to an asymptotic value φ 0 at large distances from the Sun.