Why matter effects matter for JUNO

In this paper we focus on the Earth matter effects for the solar parameter determination by the medium baseline reactor experiment such as JUNO. We derive perturbative expansions for the mixing angles θ12 and θ13 as well as the ∆m21 and ∆m31 in terms of the matter potential relevant for JUNO. These expansions, up to second order in the matter potential, while simple, allow one to calculate the electron antineutrino survival probability to a precision much better than needed for the JUNO experiment. We also quantitatively explain the shift caused by the matter effects on the solar neutrino mixing parameters ∆m21 and θ12 which do not satisfy the naive expectations and are significant given the precision that can be achieved by the JUNO experiment.


Introduction
Jiangmen Underground Neutrino Observatory (JUNO) [1] is a medium baseline reactor antineutrino experiment with ∼ 50 km kilometer baseline that is currently under construction. The primary stated goal of this experiment is to determine the neutrino mass ordering, i.e. whether the neutrino mass eigenstate with least ν e is the most massive, normal ordering or least massive, inverted ordering. One of the secondary goals of this experiment is to measure the solar neutrino mixing parameters (sin 2 θ 12 and ∆m 2 21 ) and the atmospheric ∆m 2 's with sub percent precision. In this paper we address the Wolfenstein matter effects on the measurement of the solar parameters by JUNO.
To be more precise, we consider the JUNO-like experiment, or the reactor experiment which has the same main features of JUNO in terms of its size, baseline and energy resolution with some simplifications and approximations regarding the experimental setup and/or analysis. For brevity, we refer such a experiment simply as JUNO in this paper. There is a proposal of similar experiment, RENO-50 [2], but in practice, currently, JUNO is the only reactor experiment with medium baseline which is expected to start taking data soon, within a next few years. Therefore, we focus on the JUNO experiment in this paper.
We first derive a simple perturbative expansion for the mixing angles and ∆m 2 's in matter. Using these matter mixing angles and matter ∆m 2 's one can calculate the Email addresses: amir.khan@mpi-hd.mpg.de, akhan@fnal.gov (Amir N. Khan), nunokawa@puc-rio.br (Hiroshi Nunokawa), parke@fnal.gov (Stephen J. Parke) electron antineutrino survival probability with a precision that is beyond what is needed for the JUNO experiment. The fractional precision is better than 10 −5 . The matter effects only significantly impact the measurements of the solar parameters sin 2 θ 12 and ∆m 2 21 . Previous attempts to approximate the oscillation probability relevant for experiment like JUNO including matter effects [3][4][5] seem to be either very complicated or less accurate than the ones shown in this paper. The matter mixing angles and mass-squared differences derived here are simple expansions in the matter potential up to second order.
Due to small matter effect, the effective mixing angle and mass squared difference in matter are slightly shifted from the ones in vacuum, if they are obtained by ignoring intentionally the matter effect. Certainly, the matter effect will be included when JUNO will analyze its reactor neutrino data but it is an interesting exercise to see explicitly the impact of matter effect on the solar parameter determination as studied in [5]. While it would not be possible to identity or establish the presence of the matter effect by JUNO alone, one can perform a consistency check.
In one of the previous studies, in [5], using a full χ 2 analysis, it was shown that the fractional shift for the solar parameters due to the Earth matter effect (assuming the constant matter density 2.6 g/cm 3 ) are δ(sin 2 θ 12 ) sin 2 θ 12 , δ(∆m 2 21 ) ∆m 2 21 χ 2 analysis On the other hand, by considering only the first order in matter effect, from eqs. (11) and (17) Comparing these results it is clear that the naive expectation can explain either the shift for sin 2 θ 12 or ∆m 2 21 but not both at the same time using a fixed neutrino energy. We have checked that including also higher order matter corrections does not help this situation. In this paper we explain in a semi-analytical way this apparent discrepancy we observe and confirm that our predictions are consistent with the results of χ 2 analysis of [5].
In section 2, we give the expansion for the mixing angles, sin 2 θ's and ∆m 2 's as expansions in the matter potential which are used to calculate the electron anti-neutrino survival probability with sufficient precision for the JUNO experiment. In section 3, we discuss the event rates both in vacuum and in matter taking into account the smearing for the reconstructed neutrino energy due to finite energy resolution. In section 4 we give the naive estimate of the shift in sin 2 θ 12 and ∆m 2 21 as well as a semi-analytic estimate of the shift using the reactor spectral information as well as the energy dependence of the inverse beta decay cross section. Section 5 provides the summary and conclusions. In Appendix we provide some details for the calculations of the shift given in the section 4.
Here we will use the approximate expressions derived by Denton, Minakata and Parke (DMP) [7] which depend on two quantities, first the Wolfenstein matter potential times neutrino energy, where the +(-) sign refers to the case where neutrino (antineutrino) channel is considered, G F is the Fermi constant, N e is the electron number density, and second, ∆m 2 ee ≡ cos 2 θ 12 ∆m 2 31 + sin 2 θ 12 ∆m 2 32 .
Since we are interested in reactor neutrinos to be observed by JUNO, in this work we will mainly consider antineutrino channel (considering the negative sign in eq. (5)) unless otherwise stated. First for θ 13 , we have [7] cos 2 θ 13 ∆m 2 ee cos 2θ 13 − a ∆ m 2 ee , where ∆ m 2 ee ≡ ∆m 2 ee (cos 2θ 13 − a/∆m 2 ee ) 2 + sin 2 2θ 13 , (8) and then θ 12 and ∆ m 2 21 cos 2 θ 12 ∆m 2 21 cos 2θ 12 − a ∆ m 2 21 , where ∆ m 2 21 ∆m 2 21 × (cos 2θ 12 − a /∆m 2 21 ) 2 + cos 2 ( θ 13 − θ 13 ) sin 2 2θ 12 , (10) where the effective matter potential for 1-2 sector, a , is given by and and finally ∆ m 2 31 and ∆ m 2 32 ; Performing a Taylor series expansion in the parameters a/∆m 2 ee and a/∆m 2 21 on eqs. (7)-(10) one can calculate all the oscillation parameters up to the required accuracy level. In the following, we write down the mass-squared differences and the mixing angles in the parameters up to the order accurate enough for JUNO:  Note that the solar parameters, sin 2 θ 12 and ∆ m 2 21 , are given to second order in a/∆m 2 21 whereas the atmospheric parameters sin 2 θ 13 and ∆ m 2 ee are given to only first order in a/∆m 2 ee . However, the solar parameters use ∆m 2 21 , whereas the atmospheric parameters use ∆m 2 ee . Numerically, , for sin 2 θ 13 given in Table 1, so that at E ∼ 11 MeV which is just beyond the highest energy of reactor neutrinos.
In Fig. 1 we show the matter effects on each parameter relevant for JUNO using the mixing parameters shown in Table 1. In left panels we show the ratio of the mixing parameters in matter and the corresponding ones in vacuum as a function of neutrino energy whereas in the right panels we show the fractional precision given by our approximation in eqs. (15)-(20). We note that positive (negative) energy corresponds to the case of antineutrino (neu-trino) channel since as long as the matter effect is concerned, changing the sign of neutrino energy is equivalent to changing that of the matter potential, see eq. (5).
Clearly the matter effect on solar parameters is more than an order of magnitude larger than for the atmospheric parameters. The fractional precision figure shows that the fractional difference between our approximated expressions of eqs. (15)-(20) and the exact values are less than 5 × 10 −6 for E below 10 MeV.
In Fig. 2 we show the fractiona difference between survival probability in matter and that in vacuum versus the antineutrino energy. The figure clearly shows that the most important region is between 2 to 5 MeV and the peak occurs around 3 MeV. This regime is also the most important for the precision measurement of the solar oscillation parameters.
The accuracy of our approximate probability using eq. (4) with effective mixing parameters in eqs. (15)-(20) for the energy range relevant for the reactor spectrum and JUNO baseline is shown in Fig. 3. We took the absolute difference of our approximate probability with the exact probabilities using the exact oscillation probabilities of ref. [6]. For comparison, we also show the difference between the full DMP probability [7] against the exact probability in blue color.
In addition, we show by the solid purple curve the same quantity but computed using only the first order in matter effect. We note that roughly speaking, if we consider only the statistical error, JUNO is expected to measureν e survival probabilities with an error of ∼ O(0.1) % which is just a factor of few times larger than the errors of the probability with the first order in matter effect around 4 MeV.

Event Rates
For the purpose of this paper, we consider the JUNO detector with the exposure corresponding to (36 × 20 × 5) GW·kton·years of reactor power times target mass for five years of running time. We compute the differential event rate in terms of the reconstructed neutrino energy, E rec , as where N p is the total number of target protons in the liquid scintillation detector of JUNO, L is the baseline, T is the total run time, dφ/dE is the reactor flux we take from ref. [9], P ee (L, E) is the oscillation probability, σ IBD is the total cross-section of inverse beta decay at the detector taken from ref. [10]. For simplicity and as a good approximation we ignore the small variation (of the order of ∼ 0.5%) of the distances from the JUNO detector to the 10 reactor cores and set L = 52. 5    In the left panels we show the ratios of all the mixing parameters which are relevant for medium baseline reactor neutrino experiment, in matter and vacuum, or "exact matter/vacuum", computed using the exact formulas in matter [6] as a function of neutrino energy. In the right panels we show the fractional error between exact and approximated values of mixing parameters in matter, or "| approximated matter -exact matter | /exact matter", as a function of neutrino energy where eqs. (15)-(20) based on DMP results were used for the approximated ones. We note that the positive (negative) energy corresponds to the antineutrino (neutrino) channel since as long as the matter effect is concerned, changing the sign of energy is equivalent to changing the sign of the potential in eq. (5). The difference between exact and either of the approximations would be invisible on the left panels.
are interested to estimate only the matter induced shift of the solar parameters (not the accuracy of the measurement), we believe that taking into account the variation of baselines have essentially no impact on our results. This is expected to be true because we treat both the input and fit exactly in the same way (by applying the same approximations to both input and fit). In eq. (24) G(E − E rec , δE) is the normalized Gaussian smearing function which takes into account the photon energy smearing of the detector. We define this function to be, where δE is the energy detector resolution, defined as [1] δE following the analysis done in [11]. Here, we have used that m n − m p − m e = 0.8 MeV, so that the visible energy deposited in the detector is related to the reconstructed neutrino energy as E vis ∼ E rec − 0.8 MeV which at threshold is ∼ 1 MeV. We do not consider neither the proton recoil in the IBD reaction nor the so called "non-stochastic" terms in eq. (26) (see eq. (2.11) and related description in [1]) due to the same reason mentioned before. As long as we are interested only in the shift due to matter potential, we believe that these details are not important. For our purpose, we are satisfied by the fact that by our sim-   [7]) and the exact result. We also show by the solid purple curve the same quantity but computed using only the first order in matter effect.
plified χ 2 analysis we could reproduce rather well the results shown in Fig. 4 of [5].
In the upper panel of Fig. 4 we show the event number distribution. On the lower panel, we show the fractional difference of the event rates between matter and vacuum.

Shift of solar parameters, sin 2 θ 12 and ∆m 2 21 , due to the matter effect at JUNO
As an application of the analytic formulas shown in section 2 here we try to estimate quantitatively the size of the shift due to matter effects for θ 12 and ∆m 2 21 for JUNO. Namely, we will see how much the results differ between the 2 cases when matter effect is taken into account or ignored, by considering more details than the discussions done in [5] regarding the shift. Considering JUNO's expected sensitivity, the matter effect should definitely be taken into account.
We first try to estimate the shift due to the matter effect by using analytic probability formulas considering only the first order in matter effect to compute the impact of matter effect for a given neutrino energy. We will show that by choosing appropriately the neutrino energy (∼ 3 MeV) we can predict, without performing a fit to the data, rather accurately the size of the shift due to the matter effect for the mixing angle sin 2 θ 12 but we tend to overestimate somewhat the shift for ∆m 2 21 . Alternatively, if we choose ∼2 MeV which is at the very low end of the observed neutrino spectrum, we can predict the shift in ∆m 2 21 but underestimate the shift in sin 2 θ 12 . As mentioned in the introduction, including also higher order matter corrections does not help this situation.
Later in this section, we will resolve this apparent discrepancy in a relatively simple manner, by accurately computing or predicting the shift due to the matter effect for sin 2 θ 12 and ∆m 2 21 , without performing a full χ 2 analysis. Thus, confirming the full χ 2 analysis in an independent way as long as the shift is concerned.

Naive estimation of the shift due to matter effect
First let us try to estimate analytically the expected shift due to matter effect in the determination of the solar parameters sin 2 θ 12 and ∆m 2 21 by JUNO. For the sake of illustration and simplicity, we consider only the first order in matter effect in eqs. (15) and (16), and adopt the ap-proximation that the rapid oscillation driven by the atmospheric mass squared difference is simply averaged out, since as long as the estimation of the shift due to matter effect for solar parameters is concerned, this rapid oscillation is expected to be unimportant.
Under this approximation, from eq. (4) we can obtain the survival probability with matter effect which is essentially reduced to that for 2 flavor system with some small corrections due to non-zero θ 13 , which can be expressed asP ee (E) 1 − c 4 13 sin 2 2θ 12 sin 2 ∆m 2 21 L 4E − 1 2 sin 2 2θ 13 , (27) whereθ 12 and ∆m 2 21 are given, respectively, through the following expressions, sin 2θ 12 sin 2 θ 12 1 − 2α(ρ)E cos 2 θ 12 cos 2θ 12 , or equivalently, and ∆ m 2 which are obtained by taking into account only the first order in matter effect in eqs. (15) and (16). We note that shift indicated in eqs. (29) and (30) agree, respectively, with the ones shown in eqs. (11) and (17) of [5] apart from the correction factor c 2 13 in α(ρ) due to non-zero θ 13 . From eqs. (28) and (30) we can try to estimate the expected shift due to the matter effect for JUNO by choosing some representative value of reactor neutrino energy. It seems reasonable to consider the neutrino energy corresponding to the first oscillation minimum driven by ∆m 2 21 for the JUNO baseline as follows, for L = 52.5 km. By using this energy and ρ = 2.6 g/cm 3 , we expect that if we perform a fit to the data (which inevitably include matter effect) by ignoring the matter effect, the best fitted values of sin 2 θ 12 and ∆m 2 21 would be shifted as In this work sin 2 θ fit 12 and ∆m 2 12 fit imply the mixing parameters to be obtained by fitting the data in a χ 2 analysis, which can be performed with or without taking into matter effect. Since one can not switch off the matter effect in a real data, we assume that the input data to be fitted always include the matter effect whereas it can be either included or neglected in the fit. This implies that if we fit the data by using vacuum formula, we tend to underestimate (overestimate) the value of sin 2 θ 12 (∆m 2 21 ) by these percentages using E osc. min . Now let us compare this estimation with the results which can be obtained by fitting the data. As mentioned in the introduction, the naive expectation of the shift given in eq. (33) does not agree with the the results of the full χ 2 analysis done in [5] (see Fig. 5 of this reference) which implies that the shift is δ(sin 2 θ 12 ) sin 2 θ 12 , δ(∆m 2 21 ) ∆m 2 21 χ 2 analysis which we confirmed also by our χ 2 fit performed in a similar way as done in [11], giving very similar results for the fractional differences of (sin 2 θ 12 , ∆m 2 21 ) ∼ ( -1.2, 0.20)%. By comparing the values given in eqs. (33) and (36), we observe that while the shift for the mixing angle agrees quite well between the one expected by naive prediction and with that obtained by χ 2 analysis, we see that the shift for ∆m 2 21 does not agree very well (0.33% vs 0.19%). As long as ∆m 2 21 is concerned we would need to use E ∼ 1.9 MeV to get the correct shift but this looks odd as this value is too close to the energy threshold for the inverse beta decay, and moreover, with such a energy, the agreement for the shift for the mixing angle gets worse, as mentioned in section 1 (introduction). We have checked that taking into account higher order correction of matter effect would not help to make the agreement better.
In the next subsection, we try to compute more accurately the shift of solar mixing parameters, in particular, the one for ∆m 2 21 , due to matter effect without performing a fit to the data through a detailed χ 2 analysis.

Computation of the expected shift due to matter effect
Based on the previous works such as [5] we can safely assume that in the presence of the (small) matter effect, if we try to fit the data by ignoring the matter effect, the quality of the fit would not be aggravated (compared to the case where the matter effect is correctly taken into account) but the best fit values of the mixing parameters are simply shifted from the true (correct) values by some small constants 1 .
We try to parameterize such constants by dimensionless parameters x and y as follows, sin 2 2θ fit 12 = sin 2 2θ 12 (1 + x) where x (y) is expected to take some negative (positive) value which implies that mixing angle (mass squared difference) is tend to be underestimated (overestimated), in agreement with what has been obtained previously [5].
Here we consider the shift in terms of sin 2 2θ 12 , to make the computation (shown in Appendix) simpler and then convert it in terms of sin 2 θ 12 for comparison.
Note that x and y are parameters which do not depend explicitly on neutrino energy, so that the relations in eq. (37) are different from (or not equivalent to) that given in eqs. (29) and (30). If the energy spectrum were monochromatic (i.e, E = E 0 ), we have, x = −2αE 0 and y = αE 0 , so one naively expects that x ≈ −2y also for the case of the non-monochromatic spectrum but this is turned out to be not true due to the non-trivial energy dependence in the event number distribution as we will see below.
As the derivation will be shown in Appendix, as a good approximation, the shift parameters x and y are computed as follows, where parameters a, b, c, d, f and g are computed by taking into account the normalized event number distribution in the absence of oscillation, λ(E), defined as where S(E) is the reactor neutrino spectrum taken from [9] and σ IBD (E) is the inverse-beta decay cross-section taken from [10]. Explicitly, we can compute these parameters by performing the following integrals, and with where the numerical values given in the above equations are computed using the mixing parameter and the matter density given in Table 1. By using the numerical values given in above equations we find that x = −0.317% and y = 0.198%, implying δ(sin 2 θ 12 ) sin 2 θ 12 , δ(∆m 2 21 ) ∆m 2 21 this work (2x cos 2 θ 12 / cos 2θ 12 , y), which agree quite well with the shift obtained by the χ 2 analysis done in [5] and also with the χ 2 fit performed in this work. Thus, as long as the correction due matter effects is small, it is expected that the shifted fitting parameters would be given by sin 2 θ fit 12 ≈ sin 2 θ 12 1 − 0.0065 cos 2 θ 12 cos 2θ 12 ρ 2.6 g/cm 3 , (47) ∆m 2 fit 21 ≈ ∆m 2 21 1 + 0.002 ρ 2.6 g/cm 3 .
In order to see more clearly (visually) how the shift predicted by eq. (38) agrees well with the result of the χ 2 fit, in Fig. 5, we show the results of the the χ 2 fit and shift obtained by eq. (38). In this figure we show the best fit points (indicated by the filled squares with 1 σ error bars considering only statistics for simplicity) obtained by χ 2 analysis in the plane of sin 2 θ fit 12 − ∆m 2 fit 21 for the 3 cases where (a) ρ true = ρ fit , indicated by black color, (b) ρ true = 2.6 g/cm 3 = ρ fit = 0, indicated by red color and (c) ρ true = 5.2 g/cm 3 = ρ fit = 0, indicated by blue color. For the χ 2 analysis, the probability was computed by the exact solutions with matter effect including all order of corrections.
Simultaneously in the same plot, we also show predicted shifts due to the matter effect obtained by the naive formulas in eq. (33) by the dashed magenta line and by our predictions given in eqs. (47) and (48) by the dashed green line where ρ true is varied continuously, from right to left, from 0 to 6.5 g/cm 3 and 2 open green circles and magenta squares correspond to the case where ρ true = 2.6 g/cm 3 (right open circle/square) and 5.2 g/cm 3 (left open circle/square) g/cm 3 . Note that some unlikely values of ρ true (significantly different from the reference value 2.6 g/cm 3 ) were considered just for the sake of illustration.
We first observe that naive predictions for the shift obtained by eq. (33) do not agree very well with the results of χ 2 fit because the shift for ∆m 2 21 for a given matter density is always somewhat overestimated. On the other hand, the best fit values obtained by χ 2 fit (indicated by the solid squares) agree rather well with the predictions computed by eq. (38) (indicated by open green circles), especially for the reference matter density for the JUNO baseline, ρ = 2.6 g/cm 3 . We believe that the small discrepancy we see for larger matter density, namely, for the case (c) between the best fit obtained by the χ 2 analysis (indicated by the blue square) and the predicted one (indicated by an open green circle close to the blue square) mainly comes from the higher order correction of the matter effect.

Summary and Conclusions
We have given simple, perturbative expansions for the sin 2 θ's and ∆m 2 's in the matter potential which can be used to calculate the electron anti-neutrino survival probability with precision more than necessary for the JUNO experiment. We have shown that the maximum difference between the vacuum and matter oscillation probability occurs at the solar oscillation minimum, around 3 MeV for the JUNO experiment and has a magnitude of 3.5% (4.0%) including (not including) energy resolution smearing.
Then we compare the naive matter shift for sin 2 θ 12 and ∆m 2 21 with what is obtained from a full χ 2 analysis. We have explained the apparent discrepancy between these two estimates using a semi-analytic approaches that takes into account the reactor neutrino spectrum and the energy dependence of the cross section. This provides independent confirmation of the the results of [5] using a totally different approach. Once JUNO will have real data, it could be interesting to perform also a fit to the data by ignoring the matter effect to obtain the solar mixing parameters and compare them to the ones obtained with the matter effect.  Figure 5: We show the best fit points (indicated by the filled squares with bars indicating 1 σ statistical errors only) obtained by χ 2 analysis in the plane of sin 2 θ fit 12 − ∆m 2 fit 21 for the 3 cases where (a) ρ true = ρ fit (= 0 or 2.6 g/cm 3 ), indicated by black color, (b) ρ true = 2.6 g/cm 3 = ρ fit = 0, indicated by red color and (c) ρ true = 5.2 g/cm 3 = ρ fit = 0, indicated by blue color. Simultaneously in the same plot, we also show predicted shifts due to the matter effect obtained by naive expectation in eq. (33) by the dashed magenta line and by our predictions in eqs. (47) and (48), by the dashed green line where ρ true is varied continuously, from right to left, from 0 to 6.5 g/cm 3 . Two open green circles and magenta squares correspond to the case where ρ true = 2.6 (right circle/square) and 5.2 (left circle/square) g/cm 3  where δ(sin 2 2θ 12 ) = −2αE and δ(∆m 2 21 ) = αE as, respectively, given in eqs. (29) and (30) in the presence of the matter effect.
Then the expected event number distribution (to be observed) in the presence of matter effects for small value of α(ρ)E, where α(ρ) is given in eq. (31), can be written as, N obs (E) ≈ N tot 0 λ(E)[P vac ee (E) + δP ee (−2αE, αE)], (52) where ∆ 21 ≡ ∆m 2 21 L/4E, N tot 0 is the expected total number of event in the absence of oscillation, λ(E) is the nor-malized event number distribution in the absence of oscillation given in eq. (39). Now let us assume that this event number distribution can be fitted (mimicked) by using the vacuum oscillation probability but with slightly shifted mixing parameters given in eq. (37), namely, sin 2 2θ 12 (1 + x) and ∆m 2 21 (1 + y), which do not depend explicitly on neutrino energy, as if the matter effect were absent as, where δP ee (x, y) has exactly the same functional form as given in eq. (51) but just replacing δ(sin 2 2θ 12 ) and δ(∆m 2 21 ) by x and y, respectively. Since x and y do not depend explicitly on neutrino energy, we stress that N fit (E) can not be precisely equal to N obs (E).
The dimensionless shift parameters x and y can be obtained by minimizing the following χ 2 function, with respect to these parameters, where the sum over the discrete energy bins was replaced by the integral dE due to large statistics, as done in [3], and all the systematic parameters were ignored as they are not expected to be relevant as long as the computation of the shift is concerned. After plugging N obs (E) and N fit (E) from eqs. (52) and (53) in eq. (54), and ignoring δP ee term in the denominator of the integrand of eq. (54), one can obtain The minimization condition is given by which lead to linear equations for 2 unknowns x and y, a b c d x where a, b, c, d, f and g are given in eqs. (40)-(44), and can be easily solved, as given in eq. (38). Note that the expected shift does not depend on the statistics (or the total number of events) of the experiment. In principle, taking into account higher order correction of matter effects is possible though the expressions for the solutions become more complicated.