Uncertainty in the Reactor Neutrino Spectrum and Mass Hierarchy Determination

TL;DR

This paper investigates how uncertainties in reactor neutrino spectra and fuel composition affect the determination of neutrino mass hierarchy, proposing methods to mitigate these issues with near detectors.

Contribution

It introduces a model-independent approach to handle spectrum uncertainties and analyzes the impact of fuel composition evolution on hierarchy sensitivity.

Findings

Spectrum uncertainty impacts hierarchy determination.

Large near detectors can minimize spectrum difference effects.

Fuel composition evolution can be reconstructed from near detector data.

Abstract

One of the challenges that must be overcome in order to determine the neutrino mass hierarchy using reactor neutrinos is the theoretical uncertainty in the unoscillated reactor neutrino spectrum: this is one of the reasons why, recently, it was proposed to add a near detector to the JUNO experiment. A model-independent treatment of the spectrum uncertainty will be discussed, as well as the effect that it will have on the final result. Moreover, since the neutrino spectrum depends on the chemical composition of the fuel, the spectra at the near and far detectors will be different, because they will receive neutrinos from different cores. Taking into account the time evolution of the chemical composition of the fuel in the reactor core, it is possible to reconstruct the far detector spectrum from the near detector data. We will show that the method used to reconstruct the spectrum can affect sensitivity to the mass hierarchy, however if the near detector is large enough the difference will be negligible.

Tables

Table 1. Table 1: Δ ​ χ 2 Δ superscript 𝜒 2 \Delta\chi^{2} computed with and without the approximation ( G ​ P ) ​ ( G ​ n ) 𝐺 𝑃 𝐺 𝑛 (GP)(Gn)

Asimov$\backslash$Fit	$\left(GPn\right)$	$\left(GP\right)\left(Gn\right)$
$\left(GPn\right)$	12.60	13.77
$\left(GP\right)\left(Gn\right)$	-	12.53

Equations

$$\Delta {\chi }^2={\chi }_{IH}^2-{\chi }_{NH}^2$$

$${\chi }_{MH}^2=\underset{\vec{\alpha }}{\min }\left({\chi }_{MH,F}^2(\vec{\alpha })+{\chi }_{MH,N}^2(\vec{\alpha })+P.T.(\vec{\alpha })\right)$$

$$\Delta m_{21}^2=7.53\times 10^{-5}{\text{eV}}^2\Delta m_{32,NH}^2=2.51\times 10^{-3}{\text{eV}}^2$$

$$\Delta m_{32,IH}^2=-2.56\times 10^{-3}{\text{eV}}^2{\sigma }_m=0.05\times 10^{-3}{\text{eV}}^2$$

$${\text{sin}}^{2}2{\theta }_{23}=1{\text{sin}}^{2}2{\theta }_{13}=0.083{\text{sin}}^{2}2{\theta }_{12}=0.851$$

$${\chi }_{MH,D}^2(\vec{\alpha })=\sum _i\frac{(N_{exp,D,i}-N_{Th,D,i}(MH,\Delta m_{32}^2,{\beta }_i){)}^2}{N_{exp,D,i}}+P.T.$$

$$P.T.=\frac{(\Delta m_{32}^2-\Delta m_{32,IH}^2{)}^2}{{\sigma }_m^2}$$

$$N_{MH,F,i}={\int }_{E_i}^{E_{i+1}}\text{d}E\int \text{d}{E}^{\prime }G(E;E^{\prime },{\sigma }^2(E^{\prime }))P_{MH}(E^{\prime },L)\sum _K{\varphi }_K(E^{\prime })$$

$$N_{MH,F,i}=\sum _j{G}_{ij}{P}_{MH,j}^{\prime }\sum _K{n}_{K,j}^{\prime }$$

$$G_{ij}=G(E_i+\Delta E/2;E_j+\Delta E/2,{\sigma }^2(E_j+\Delta E/2))\Delta E$$

$$P_{MH,j}^{\prime }=P_{MH}(E_j+\Delta E/2,L)n_{K,j}^{\prime }={\int }_{E_j}^{E_{j+1}}\text{d}E{\varphi }_K(E)$$

$$N_{N,i}=\frac{1}{\mathcal{R}}\sum _j{G}_{ij}{n}_{K,j}^{\prime }$$

$$n_{R,i}^{\prime }\propto n_{K,i}^{\prime }\forall R,K$$

$$N_{F,i}$$

$$N_{N,i}$$

$$\mathcal{R}=\frac{L_N^2}{L^2}\frac{M_F{W}_F}{M_N{W}_N}$$

$$\Delta {\chi }^2={\chi }_{IH}^2=\underset{{\beta }_i^{\prime },\Delta m_{32}^2}{\min }\left(\sum _i\frac{({\sum }_j{G}_{ij}\Delta P_j^{\prime }{n}_j^{\prime }-G_{ij}{P}_{IH,j}^{\prime }{\beta }_j^{\prime }{)}^2}{{\sum }_j{G}_{ij}\Delta P_{NH,j}^{\prime }{n}_j^{\prime }}+\sum _i\frac{({\sum }_j{G}_{ij}{\beta }_j^{\prime }{)}^2}{\mathcal{R}{\sum }_j{G}_{ij}{n}_j^{\prime }}+P.T.\right)$$

$$\sum _j{G}_{ij}{P}_{MH,j}^{\prime }{n}_j^{\prime }=(\sum _j{G}_{ij}{P}_{MH,j}^{\prime })(\sum _j{G}_{il}{n}_l^{\prime })$$

$$P_i=\sum _j{G}_{ij}{P}_j^{\prime }{n}_i=\sum _j{G}_{ij}{n}_j^{\prime }{\beta }_i=\sum _j{G}_{ij}{\beta }_j^{\prime }$$

$$\Delta {\chi }^2=\underset{{\beta }_i,\Delta m_{32}^2}{\min }\left(\sum _i\frac{(\Delta P_i{n}_i-P_{IH,i}{\beta }_i{)}^2}{P_{NH,i}{n}_i}+\sum _i\frac{{\beta }_i^2}{\mathcal{R}{n}_i}+P.T.\right)$$

$$\Delta {\chi }^2={\chi }_{IH}^2=\underset{\Delta m_{32}^2}{\min }\left(\sum _j\frac{(\Delta P_i{n}_i{)}^2}{P_{NH,i}{n}_i+\mathcal{R}{P}_{IH,i}^2{n}_i}+P.T.\right)$$

$$(G{P}^{\prime }{n}^{\prime })\to (G{P}^{\prime })(G{n}^{\prime })G{P}^{\prime }\to PG{n}^{\prime }\to nG{\beta }^{\prime }\to \beta$$

$$\sum _K{n}_{K,i}=\rho n_{Pu,i}+(1-\rho )n_{Ur,i}$$

$$N_{F,i}=P_{IH,i}\left(\rho (n_{Pu,i}+{\beta }_{Pu,i})+(1-\rho )(n_{Ur,i}+{\beta }_{Ur,i})\right)$$

$$n_{N,\gamma ,i}=\frac{1}{\mathcal{R}}\left({\rho }_{\gamma }(n_{Pu,i}+{\beta }_{Pu,i})+(1-{\rho }_{\gamma })(n_{Ur,i}+{\beta }_{Ur,i})\right)$$

$$\mathcal{R}=N_t\frac{L_N^2}{L^2}\frac{M_F{W}_F}{M_N{W}_N}$$

$$N_{F,i}=P_{IH,i}\mathcal{R}\sum _{\gamma }{\alpha }_{\gamma }{n}_{N,\gamma ,i}=\left(n_i+\rho {\beta }_{Pu,i}+(1-\rho ){\beta }_{Ur,i}\right)$$

$$\sum _{\gamma }{\alpha }_{\gamma }=1\sum _{\gamma }{\alpha }_{\gamma }{\rho }_{\gamma }=\rho$$

$$\Delta {\chi }^2=$$

$$\underset{\vec{\beta },\Delta m_{32}^2}{\min }\sum _i\left(\frac{(\Delta P_i{n}_i-P_{IH,i}(\rho {\beta }_{Pu,i}+(1-\rho ){\beta }_{Ur,i}){)}^2}{P_{NH,i}{n}_i}+\sum _{\gamma }\frac{({\rho }_{\gamma }{\beta }_{Pu,i}+(1-{\rho }_{\gamma }){\beta }_{Ur,i}{)}^2}{{\mathcal{R}}^2{n}_{N,\gamma ,i}}+P.T.\right)$$

$${\beta }_{Ur,i}\to {\tilde{\beta }}_{Ur,i}=(1-\rho )P_{IH,i}{\beta }_{Ur,i}{\beta }_{Pu,i}\to {\tilde{\beta }}_{Pu,i}=\rho P_{IH,i}{\beta }_{Pu,i}$$

$$\Delta {\chi }^2=\underset{\vec{\beta },\Delta m_{32}^2}{\min }\sum _i\left(\frac{(\Delta P_i{n}_i-{\tilde{\beta }}_{Pu,i}-{\tilde{\beta }}_{Ur,i}){)}^2}{{\sigma }_{Exp,i}^2}+C_{Pu,i}{\tilde{\beta }}_{Pu,i}^2+C_{Ur,i}{\tilde{\beta }}_{Ur,i}^2+2C_{Mix,i}{\tilde{\beta }}_{Ur,i}{\tilde{\beta }}_{Pu,i}+P.T.\right)$$

$$C_{Pu,i}=\sum _{\gamma }\frac{{\rho }_{\gamma }^2}{{\rho }^2{\mathcal{R}}^2{P}_{IH,i}^2{n}_{N,\gamma ,i}}{C}_{Ur,i}=\sum _{\gamma }\frac{(1-{\rho }_{\gamma }{)}^2}{(1-\rho {)}^2{\mathcal{R}}^2{P}_{IH,i}^2{n}_{N,\gamma ,i}}$$

Full text

1 Introduction

In the next decade several experiments will attempt to measure the neutrino mass hierarchy [1, 2, 3, 4, 5, 6]: some experiments are already taking data [7, 8], and preliminary results seems to point towards the normal hierarchy [9], however we are far from a clear determination. In reactor neutrino experiments, like JUNO [6], the mass hierarchy is determined by studying the interference between the 1-3 and 2-3 oscillations at intermediate baselines. There are several challenges that must be overcome in this kind of experiment: an excellent energy resolution is needed ($\Delta E/E\leq 3\%/\sqrt{E}$, where the energy is expressed in MeV), and the systematic errors on the energy reconstruction must be strongly constrained.

Another problem that must be faced is the uncertainty in the reactor neutrino spectrum: a few years ago, with the measurement of the “5 MeV bump”, it became clear that the theoretical models for reactor neutrinos are not in agreement with the experimental data; moreover all the experimental data currently available are obtained with an energy resolution considerably worse than the one required for the mass hierarchy determination (at Daya Bay, $\Delta E/E\simeq 7\%/\sqrt{E}$). It is possible (and there are also some hints in this sense also from ab initio calculations) that the reactor neutrino spectrum is not smooth as it is assumed to be in the current theoretical models, but there is a “fine structure”, present but currently undetected, that could affect the determination of the mass hierarchy. This is one of the reasons why, recently, it was announced that the JUNO experiment will have a near detector, called JUNO-TAO (Taishan Antineutrino Observatory) [10, 11], with mass of the order of tons. In this paper we will discuss how the uncertainty on the theoretical spectrum (constrained by the data from the near detector) can affect the mass hierarchy determination.

An additional complication arises from the fact that the JUNO far detector will receive neutrinos from two different power plants, Taishan (2 cores) and Yangjiang (6 cores), while the near detector will be able to measure the spectrum only of one of Taishan’s cores; since the reactors in the two complexes are of different model and generation (EPR, Gen. III at Taishan, CPR-1000, Gen. II at Yangjiang), there will be a difference between the unoscillated spectra at the near and far detector. We will discuss how to parametrize the spectra in this case and, assuming that the reactor neutrino spectrum depends only on the chemical composition of the fuel, we will show that from the time evolution of the near detector spectrum it is possible to reconstruct a spectrum with a generic chemical composition. Some of the results in this paper were also presented at the PhysStat-$\nu$ meeting at CERN, in 2019 [12].

2 Mass Hierarchy From Reactor Neutrinos

More than 15 years ago Petcov and Piai [13] pointed out that it is possible to determine the neutrino mass hierarchy from reactor neutrinos, by studying the oscillation probability at intermediate baselines. The main problem (which is somehow related to most if not all the challenges that must be faced in this kind of experiment) is that there is a strong degeneracy between a change of hierarchy and a shift in $\Delta m_{32}^{2}$.

Indeed, as is shown in Fig. 1, if we try to fit the normal hierarchy spectrum assuming the inverted hierarchy, it is possible to obtain the same high-energy oscillations by just shifting the best fit value for the inverted hierarchy by less than 1.5 $\sigma$’s (using the best fit values and the precision reported in [14]); it should be underlined that, since usually experiments cannot measure directly $\Delta m_{32}^{2}$ or $\Delta m_{31}^{2}$ but only an effective mass which depends on the hierarchy, from global fits we get two different values for $\Delta m_{32}^{2}$, one for each hierarchy; this makes more difficult to break this degeneracy, even if constrains from other experiments are taken into account.

The main consequences of the degeneracy can be summarized as follows:

• •

Using a $\chi^{2}$ test the high-energy part of the spectrum is needed to break the degeneracy; the hierarchy can then be determined by comparing the oscillations at low-energy

• •

However, since these oscillation are very fast, an excellent energy resolution is needed at low energies: the goal for JUNO is to achieve an energy resolution of $3\%/\sqrt{E}$; even with such resolution, the energy smearing will significantly decrease the amplitude of the oscillations. From Fig. 2 it can be read that, using the value of $\Delta m_{32}^{2}$ that would minimize the $\chi^{2}$, the difference between the spectra is of the order of $1-4\%$ and that most of the information on the hierarchy comes from the region between 2 and 4 MeV.

• •

Since the hierarchy is determined by a small shift in the position of the low-energy peaks, even a very small systematic error on the energy reconstruction could significantly affect the final result. For this reason, this kind of experiment requires also very good constraints on the non-linearity [15, 16, 17].

2.1 Fine Structure

In nuclear reactors neutrinos are not produced in the fissions, but from the beta decay of the unstable isotopes created. Even if, in theory, it is possible to calculate explicitly the neutrino spectrum from first principles (ab initio calculations), it is not yet possible to reach a good precision in this way, due to the large number of isotopes and different decays that must be taken into account (of the order of $10^{3}$ and $10^{4}$, respectively). The flux model most commonly used in reactor neutrino experiments is the Huber+Muller model [18, 19], where the neutrino spectrum is obtained using the conversion method, starting from the beta spectrum of the isotopes created in the fission process. Its predictions are however in tension with the measurement of the “5 MeV bump.”

Even if it’s not possible to calculate precisely the reactor neutrino spectrum from ab initio calculations, it is nonetheless possible to obtain some important information from this approach. Due to the Coulomb enhancement of low-energy electrons in beta decay, there is a sharp cut-off at the end-point of the energy spectrum of neutrinos from beta decay. Since, in nuclear reactors, the total neutrino spectrum is given by the sum of a very large number different beta decays, this leads to the presence of sawtooth-like features in the spectrum (also called “fine structure”) [20, 21]. The fact that these have not yet been observed in reactor neutrino experiments is not surprising, since all the experimental data currently available are obtained with an energy resolution considerably lower than what required for the mass hierarchy determination. Indeed, the fine structure is expected to be very small, of the order of 1-2%, and such small fluctuations would be completely suppressed by the finite energy resolution at past reactor neutrino experiments.

However, due to the degeneracy shown in Fig. 1, even a small uncertainty on the spectrum can affect the sensitivity to the mass hierarchy [22]. It was argued in [23, 24] that, using the Fourier transform analysis, the fine structure will have very little effect on the mass hierarchy determination. The reason for this is that the Fourier transform is sensitive only to the global oscillation behavior: the fine structure, in order to affect the determination of the hierarchy, should have a frequency similar to the one of 1-3 and 2-3 oscillations, which the authors claim is unlikely, since the two processes are based on completely different physical mechanisms. Clearly the impact of each systematic error depends on the analysis method and so even if the fine structure has less effect on the Fourier transform analysis this does not imply that it can be neglected using, for example, a $\chi^{2}$ analysis.

It is also important to notice that there are other serious problems related to the determination of the neutrino mass hierarchy using the Fourier transform. As mentioned before, one of the most difficult challenges in this kind of experiment is to control the systematic errors due to the non-linearity: with a $\chi^{2}$ analysis, usually a precision of at least $1\%$ is required (which is not trivial to achieve in such a large detector); however using the Fourier transform an even better precision should be required [16]; in [23] it is mentioned that for this kind of analysis the precision on the unknown non-linearity should be at least $0.5\%$. Moreover the Fourier transform is more sensitive to the high-energy part of the spectrum, where the flux is very low and the background is more relevant. Using a $\chi^{2}$ analysis it is possible to simply ignore the high-energy events, since they carry very little information on the mass hierarchy, however if we use the Fourier analysis, a high-energy cut will introduce a spurious and unwanted dependence on $\Delta m_{32}^{2}$, as was shown in [25], introducing a new source of error. These are some of the reasons why, in the last years, this method was used very rarely in the literature, and in almost all the papers on the topic a $\chi^{2}$ analysis was used.

3 $\chi^{2}$ Analysis

Our test statistic is defined as:

[TABLE]

Here $\chi^{2}_{IH(NH)}$ is the $\chi^{2}$ statistic where the inverted (normal) hierarchy is assumed. We can write

[TABLE]

where $MH=NH,IH$, $\vec{\alpha}$ is a vector that contains all the pull parameters (if present) that must be minimized, $\chi^{2}_{MH,F(N)}$ is the contribution to the $\chi^{2}$ of the far (near) detector and $PT(\vec{\alpha})$ represents the eventual penalty terms related to the pull parameters. In this paper we will introduce one pull parameter for each energy bin, called $\beta_{i}$, to parametrize the uncertainty on the reactor neutrino spectrum: there will be no correlation between two of these parameters from different energy bins however, for a given bin, the pull parameter used at the near and far detector will be the same. Only one other pull parameter will be used, to take into account the uncertainty on the value of $\Delta m_{32}^{2}$, since the uncertainties over the other mixing parameters, the total flux renormalization, etc… have very little effect on the final result. Even if we introduce one pull parameter for each energy bin, the number of degrees of freedom is still positive, because now we have two data points per bin, one for each detector.

For $\beta_{i}$ no penalty term will be introduced. This means that present theoretical and experimental knowledge of the spectra is not used, as existing constraints are negligible compared with those that will come from the near detector’s data. $\sigma_{m}$ will be $0.05\times 10^{-3}$ eV2, however it should be remembered that the best fit values for $\Delta m_{32}$ are different if the normal or inverted hierarchy is assumed, so this will also affect the penalty term. The mass hierarchy is assumed to be normal, we use the following values for the mixing parameters (all the values are taken from the PDG [14], except for $\textrm{sin}^{2}2\theta_{23}$, which is assumed to be 1):

[TABLE]

For the far detector, we will consider an effective mass of 20 ktons, while for the near detector, if not otherwise specified, we consider 3 tons, similar to the mass of JUNO-TAO as reported in [10].

It should be noted that, even if we choose parameters similar to the ones of the JUNO experiment, the focus of this study is to investigate the effect of uncertainty in the reactor neutrino spectrum on the mass hierarchy determination, not to discuss the precision that can be achieved in a specific experiment. For this reason, in order to simplify the calculations, we do not take into account effects that are not related with this problem, even if they can strongly affect the expected $\Delta\chi^{2}$; for example, we assume that all the cores are at the same distance from the far detector, while in reality this is not the case (and this yields a sizeable correction to the final result [26]). Moreover, we do not consider the systematic error related to the non-linearity, which will also strongly affect the sensitivity of the experiment. All these effects are not related to the topic discussed here, so they will be ignored. As a result the $\Delta\chi^{2}$ reported are significantly overestimated and should not be taken as representative of the precision that can be achieved at JUNO.

The expression for $\chi^{2}_{MH,D}$ (where $D=N,F$) is:

[TABLE]

where the index $i$ indicates the energy bin, and runs from 1 to $N_{E}$ and $P.T.$ is the penalty term for $\Delta m_{32}^{2}$, namely

[TABLE]

In this paper we will use the Asimov data set, where the expected value of $\Delta\chi^{2}$ for a given experiment is computed using the expected number of events instead of $N_{exp,D,i}$. As a consequence, $\chi^{2}_{NH}=0$ (since we assumed the hierarchy to be normal) and $\Delta\chi^{2}=\chi^{2}_{IH}$. For this reason, from here on we will drop the subscripts $exp$ and $Th$: the two sets can be recognized because the Asimov data set will be calculated assuming the normal hierarchy, the fit assuming the inverted. It should be noticed that, in principle, one could use two different spectra for the fit and for the Asimov data set, because we cannot assume the model we use is the correct one, however since there is a pull parameter for every bin, the difference can be eliminated just by a shift of $\beta_{i}$; which means that we can safely use the same unoscillated spectrum for both datasets.

The expected number of events at the far detector is given by

[TABLE]

where $K$ indicates the different reactor cores and it is summed over all of them, $E$ is the visible energy, $E^{\prime}$ the real energy of the neutrino, $G(E;E^{\prime},\sigma^{2}(E^{\prime}))$ is a Gaussian distribution with mean $E^{\prime}$ and standard deviation $\sigma(E^{\prime})$ that describe the finite energy resolution of the detector, $\sigma(E^{\prime})=3/100\sqrt{E^{\prime}}$, $L$ is the baseline, $P(E^{\prime},L)$ is the oscillation probability $\phi_{K}(E^{\prime})$ the unoscillated neutrino spectrum (taking into account the emission spectrum, the chemical composition of the fuel at the core, the cross section, the detector mass and the eventual detector efficiency; in principle $\phi_{K}$ depends also from the baseline, due to the geometrical factor, however this dependence is left implicit). Here we assumed that all the cores are at the same distance from the detector; this is not necessarily true, however the generalization is quite straightforward, and it will be briefly discussed in the next section.

Since, due to the degeneracy with $\Delta m_{32}^{2}$, the mass hierarchy must be determined by studying relatively small differences between very fast oscillations, the dimension of the energy bins must be small, otherwise these differences would be averaged out by the size of the bin.

Fig. 3 shows the expected value of $\Delta\chi^{2}$ as a function of the number of bins (considering the spectrum in the region between 1.5 and 8.5 MeV); in the following calculations we will use 10 keV energy bins (i.e. we will consider 700 bins in the said interval). For this reason (3.5) can be computed using an discrete approximation, as a product of matrices and vectors; assuming that $E_{i+1}-E_{i}=\Delta E$ for every $i$,

[TABLE]

where here and throughout the paper lower-case Roman indices indicate the energy bin, while upper-case Roman indices indicate the different reactor cores, and

[TABLE]

At the near detector we can use two approximations to simplify the calculations: first of all we can drop $P^{\prime}_{MH,j}$, since the neutrino oscillation can be safely neglected, second the detector will receive neutrinos only from one core. We can write

[TABLE]

where $\mathcal{R}$ is a renormalization factor that takes into account the differences in mass, neutrino flux and baseline between the near and far detector. We also assume that the far and near detectors have the same energy resolution. Notice that here $K$ is not summed but indicates a specific reactor core.

3.1 Example: Same chemical composition

Now we will consider a simple example, where we assume that all the cores have the same chemical composition, hence

[TABLE]

Here we have $\propto$ and not $=$ because the reactors can still have different thermal power.

We can express the expected number of events (writing explicitly the dependence on the pull parameters) as

[TABLE]

Here $\mathcal{R}$ is defined as

[TABLE]

where $L$ and $L_{N}$ are the baselines of the far and near detector, respectively, $M_{N(F)}$ and $W_{N(F)}$ are the mass of the near (far) detector and the thermal power of the core(s) seen by each detector. Taking into account (3.3), we can rewrite (3.2) as

[TABLE]

With perfect energy resolution (i.e. if $G_{ij}=\delta_{ij}$), the minimization over $\beta^{\prime}_{i}$ would be trivial, because the pull parameters in each bin can be minimized separately. However the finite energy resolution mixes all the pull parameters from different energy bins, and in order to obtain an analytical solution it is necessary to compute the inverse of $G_{ij}$, which would lead to huge numerical errors. It is still possible to find the minimum numerically, however due to the large number of parameters it would be quite time-consuming and it would require high computational power.

It is possible to obtain analytical results by using the following approximation:

[TABLE]

The exact results will be sometimes indicated as $(GP^{\prime}n^{\prime})$ and the approximation as $(GP^{\prime})(Gn^{\prime})$. With this approximation we are performing the convolution with the Gaussian separately for the spectrum and the oscillation probability. Assuming that the energy resolution of the near and far detector is the same, this is equivalent to considering the spectrum seen at the near detector as the ”real” unoscillated spectrum, and the expected spectrum at the far detector is obtained by multiplying the near detector spectrum with the oscillation probability convoluted with a Gaussian: what is lost in this picture is that, following the exact procedure, in the Gaussian convolution the oscillation probability should be weighted by the value of the spectrum at that particular point, however since the mass hierarchy determination requires a very good energy resolution and the Gaussian must be sharply peaked on the real energy, the error due to this approximation is quite small. Does this mean that this is a reasonable approximation? The answer is “Not really, but…”.

As can be seen in Fig. 4, except for the very low energy region (where the spectrum is sharply increasing, this region does not carry important information for the hierarchy determination), the difference between the exact calculation and the approximation is less than $0.5\%$; however this error is (albeit smaller) of the same order of magnitude as the difference between the hierarchies, this means that it could affect the final result. It is however still possible to use this approximation if it’s employed to compute both the fitting model and the Asimov data set: in this way, indeed, the $0.5\%$ error would be not on the spectrum, (i.e. $GP^{\prime}n^{\prime}$), but on the difference between the spectra ($G\Delta P^{\prime}n^{\prime}$), which is already very small and it would be a second order effect, that can be neglected. On the other hand, however, this means that this technique can be used only in numerical simulations, while in order to analyze the actual results of an experiment the minimization must be performed numerically.

In Tab. 1 we show $\Delta\chi^{2}$ computed in all of the possible cases, assuming a perfect knowledge of the unoscillated spectrum (i.e. setting $\beta_{i}=0$) and taking into account only the contribution of the far detector. We remind the reader that if the approximation is used only for the fit but not for the Asimov dataset, $\chi^{2}_{NH}$ is not zero anymore, and it must be taken into account.

We can see that if the approximation is used only in the fit, the difference with respect to the exact case is quite large ($\simeq 10\%$), while if $(GP^{\prime})(Gn^{\prime})$ is used both for the fit and for the Asimov dataset the difference is around $0.5\%$.

We can now rename

[TABLE]

The $\Delta\chi^{2}$ now reads

[TABLE]

The minimization over $\beta_{i}$ is now trivial, and the result is

[TABLE]

4 Chemical Composition

The neutrino spectrum depends on the chemical composition of the fuel. There were several works in the last few years on this topic, studying, for example, the dependence on the chemical composition of the 5 MeV bump or of the reactor neutrino anomaly [27, 28, 29]. Since at JUNO the far detector will receive neutrinos from two power plants (with reactors of different model and generation), while the near detector will see only one core, it is quite likely that the two unoscillated spectra will be different. However, the chemical composition will evolve with time, so at the near detector it will be possible to sample different compositions during the experiment: is it possible to use this information to reconstruct the spectrum at the far detector? We will consider only the contribution of two isotopes to the spectrum 239Pu and 235U, that together produce almost 90$\%$ of the neutrinos detected at Daya Bay; increasing this number would complicate the calculations but not give any additional insight on the problem. First of all we will use the same approximation employed in the previous section:

[TABLE]

We can write the unoscillated spectrum at the near detector in Eq. (3.6) as:

[TABLE]

where $\rho$ is the fraction of the fissions due to 239P, averaged over all the cores, and $n_{Pu(Ur),i}$ is the unoscillated spectrum generated by 239Pu (235U); we used “Ur” instead of the chemical symbol for Uranium in order to avoid any confusion with the indices used to indicate the reactor cores.

We must now to introduce two sets of pull parameters $\beta_{i}$, one for each isotope. Writing them explicitly we have

[TABLE]

We divide the spectrum at the near detector in $N_{t}$ time bins, all of the same size; in each time bin $\gamma$, the fraction of fissions due to 239Pu will be $\rho_{\gamma}$. In order to be able to reconstruct the unoscillated spectrum at the far detector, $N_{t}$ must be equal or larger than the number of isotopes. We have

[TABLE]

where here and in the rest of the paper the Greek indices ($\gamma$, in this case) indicate the time bins. $\mathcal{R}$ is now defined as

[TABLE]

the only difference with respect to Eq. (3.11) is the presence of the number of time bins $N_{t}$.

We can now write the unoscillated spectrum at the far detector as a linear combination of the spectra measured at the near detector in the different time bins

[TABLE]

where $n_{i}=\rho n_{Pu,i}+(1-\rho)n_{Ur,i}$ and, in the last step, we imposed the following conditions on the coefficients $\alpha_{\gamma}$:

[TABLE]

There are two important observations that must be made here: first, if $N_{t}$ is larger than the number of isotopes considered, the choice of $\alpha_{\gamma}$ is not unique; however we can notice that $\alpha_{\gamma}$ do appear not on the right side of Eq. (4.6). Indeed, since the $n_{i}$’s and $\beta_{i}$’s are multiplied by the same coefficients, the conditions (4.7) ensure that $\Delta\chi^{2}$ does not depends on our choice of $\alpha_{\gamma}$, which will not affect the final result. Second, it is not always possible to express the expected number of events at the far detector in a simple expression such as Eq. (4.6) (for example, this is not possible when the interference is taken into account), however in order to use this approach it is sufficient that we are able to reconstruct the unoscillated spectrum for each isotope as a linear combination of the near detector data.

The $\Delta\chi^{2}$ now reads

[TABLE]

where $\vec{\beta}=(\beta_{Pb,i},\beta_{Ur,i})$. Using the following transformation for the coefficients $\beta$:

[TABLE]

and expanding the the second term in the $\Delta\chi^{2}$ we can write:

[TABLE]

where

[TABLE]

It is now possible to perform analytically the minimization over the $\tilde{\beta}_{Pu,i}$’s and $\tilde{\beta}_{Ur,i}$’s, we have:

[TABLE]

To obtain these results we have assumed that the energy resolution is the same at the near and far detector; it is worth noticing that, while in the case of identical chemical composition it is possible to get an analytical result for any energy resolution, if the unoscillated spectra at near and far detector are different it is possible to write explicitly the solution only if $\sigma_{F}>\sigma_{N}$, otherwise it is necessary to compute the inverse of $G_{ij}$.

Assuming the change of the chemical composition will be constant in time, we have that $\rho_{\gamma}$ is

[TABLE]

where $\gamma=1,\dots,N_{t}$. We assume that, at the near detector, $\rho_{min}=0.25$ and $\rho_{max}=0.35$ (similar to that reported in [29], regarding the fuel evolution at Daya Bay; it is worth noticing that the data reported there involve also an average over several cores, so in the case of only one core the variation could be larger).

In Fig. 5 we show $\Delta\chi^{2}$ as a function of the number of time bins: it is possible to see that increasing $N_{t}$ will increase the precision, however with a detector in the ton range the difference is minimal. In Fig. 6 we can see the expected value of $\Delta\chi^{2}$ as a function of the different (average) chemical composition seen by the far detector and of the mass of the near detector; $N_{t}=2$.

The inclusion of the interference effect is quite straightforward: if the reactor cores are not all at the same distance from the far detector, the expected number of events can be written as

[TABLE]

Here $W_{K}$ and $L_{K}$ are the power and the baseline of the $K$-th core, $P(L_{K})_{NH,i}$ is the oscillation probability in the $i$-th energy bin, where now the dependence on the baseline is written explicitly. To be consistent with the notation, $n_{Pu,i}$ and $n_{Ur,i}$ must also be adequately normalized, i.e. they represent the unoscillated spectrum that would be measured at the far detector from a 1 GW core (with the appropriate chemical composition) 1 m away, if the power and the baseline are expressed in GW and m, respectively.

We can see that now we cannot define a generic “oscillation probability”, since it will be different for each isotope: we can collect of the factors multiplying the unoscillated spectrum by writing

[TABLE]

which, of course, do not correspond to an oscillation probability. We can follow the same procedure as before to calculate analytically the minimization over the $\beta$’s, however Eq. (4.9) now reads

[TABLE]

$\Delta\chi^{2}$ will still have the same expression as in Eq. (4.12), however now the $C$’s and $\sigma_{Exp,i}$ will be defined as

[TABLE]

4.1 Energy Reconstruction

We conclude this section with an observation: an approach similar to the one described here could be used also for a model-independent treatment of the systematic errors due to the energy non-linearity. We can indeed rewrite Eq. (3.5) as

[TABLE]

If there is an unknown systematic error in the energy reconstruction, then the visible energy is related to the real energy of the neutrino by the formula

[TABLE]

Eq. (4.18) becomes

[TABLE]

where $\epsilon_{i}=\epsilon(E_{i})$ and all the terms proportional to $\epsilon_{i}(E_{i+1}-E_{i})$ are neglected. It is now possible to use the same approach followed in the previous section to take into account the uncertainty due to the non-linearity, without the need to rely on any particular model. It is important to underline, however, that the main difficulty in this procedure is to estimate the penalty term related to $\epsilon(E)$.

To clarify this point let consider the following example: let us assume that, in the case discussed in the previous section, there is a systematic error in the unoscillated spectrum which is constant for every energy bin, i.e. $\beta_{i}=c$. Since the penalty term for $\beta_{i}$ comes from the measurement at the near detector, and the error in that measurement goes like $1/\sqrt{N_{F,i}}$, it is easy to prove that if we change the size of the bins, the sum of the penalty terms for all the $\beta$’s would remain constant. If, in the case of non-linearity, we assume naively that the uncertainty for each parameter $\epsilon_{i}$ is $1\%$ (for example), it is also easy to see that, in case of constant bias, the total penalty term would be proportional to the number of energy bins, which is clearly an undesirable feature. This means that the penalty term function could be considerably more complicated, and it could also depend on the choice of the bin size, as well as on the calibration procedure. This also introduces some sort of model-dependence in our analysis (because the final result depends on how the error bars for the non-linearity are chosen), however this is unavoidable, regardless of the method chosen; the main advantage of this approach is that it is not necessary to chose a specific model to parametrize $\epsilon(E)$.

5 Conclusion

The uncertainty in the reactor neutrino spectrum could affect the mass hierarchy determination, however we have shown that it is possible to reduce significantly the effect on the final results using the data from a near detector, with an approach that does not rely on any theoretical model; if the mass of the near detector is in the ton range, the loss of precision is negligible. Even if the chemical compositions of the reactors seen by the near and far detector are different, it is possible to reconstruct the unoscillated spectrum at the far detector by studying the time evolution of the near spectrum. One of the consequences is that the data collected at JUNO-TAO could be invaluable not only for JUNO but also for other reactor neutrino experiments, that could use those data, where the spectrum is measured with an excellent energy resolution, to extrapolate the spectrum that will be seen in their detectors, even if the chemical composition of the fuel is different.

Acknowledgement

EC thanks Bryce Littlejohn and Andrew Conant for the useful discussions and comments. JE is supported by the CAS Key Research Program of Frontier Sciences grant QYZDY-SSW-SLH006 and the NSFC MianShang grants 11875296 and 11675223. EC is supported by NSFC Grant No. 11605247. JE and EC also thank the Recruitment Program of High-end Foreign Experts for support.