Mapping reactor neutrino spectra from TAO to JUNO

The Jiangmen Underground Neutrino Observatory (JUNO) project aims at probing, at the same time, the two main frequencies of three-flavor neutrino oscillations, as well as their interference related to the mass ordering (normal or inverted), at a distance of ~53 km from two powerful reactor complexes in China, at Yangjiang and Taishan. In the latter complex, the unoscillated spectrum from one reactor core is planned to be closely monitored by the Taishan Antineutrino Observatory (TAO), expected to have better resolution (x 1/2) and higher statistics (x 30) than JUNO. In the context of neutrino energy spectra endowed with fine-structure features from summation calculations, we analyze in detail the effects of energy resolution and nucleon recoil on observable event spectra. We show that a model spectrum in TAO can be mapped into a corresponding spectrum in JUNO through appropriate convolutions. The mapping is exact in the hypothetical case without oscillations, and holds to a very good accuracy in the real case with oscillations. We then analyze the sensitivity to mass ordering of JUNO (and its precision oscillometry capabilities) assuming a single reference spectrum, as well as bundles of variant spectra, as obtained by changing nuclear input uncertainties in summation calculations from a publicly available toolkit. We show through a chi-squared analysis that variant spectra induce little reduction of the sensitivity in JUNO, especially when TAO constraints are included. Subtle aspects of the statistical analysis of variant spectra are also discussed.

The Jiangmen Underground Neutrino Observatory (JUNO) project aims at probing, at the same time, the two main frequencies of three-flavor neutrino oscillations, as well as their interference related to the mass ordering (normal or inverted), at a distance of ∼ 53 km from two powerful reactor complexes in China, at Yangjiang and Taishan. In the latter complex, the unoscillated spectrum from one reactor core is planned to be closely monitored by the Taishan Antineutrino Observatory (TAO), expected to have better resolution (×1/2) and higher statistics (×30) than JUNO. In the context of ν energy spectra endowed with fine-structure features from summation calculations, we analyze in detail the effects of energy resolution and nucleon recoil on observable event spectra. We show that a model spectrum in TAO can be mapped into a corresponding spectrum in JUNO through appropriate convolutions. The mapping is exact in the hypothetical case without oscillations, and holds to a very good accuracy in the real case with oscillations. We then analyze the sensitivity to mass ordering of JUNO (and its precision oscillometry capabilities) assuming a single reference spectrum, as well as bundles of variant spectra, as obtained by changing nuclear input uncertainties in summation calculations from a publicly available toolkit. We show through a χ 2 analysis that variant spectra induce little reduction of the sensitivity in JUNO, especially when TAO constraints are included. Subtle aspects of the statistical analysis of variant spectra are also discussed.
Another layer of complexity, pointed out in summation calculations of the neutrino energy spectrum, is the presence of sawtooth-like substructures, as expected from Coulomb effects in individual β decays [54,55]. These fine-structure features have not been observed within the resolution of current experiments, with the possible exception of a hint discussed in [56]. Observing or constraining (at least a few) prominent substructures in future, high-resolution and high-statistics reactor neutrino experiments, would be beneficial both for nuclear spectroscopy (allowing to pin down the spectral contributions of specific fission products [56]), and for neutrino oscillometry (reducing small-scale fuzzy uncertainties that might affect the JUNO sensitivity to mass ordering [57]). Although the latter benefit may be marginal if one assumes "known" substructures from nuclear databases [58,59], the observation of unexpected spectral anomalies at large energy scales (normalization and bump issues) provides a warning about the possible emergence of "unknown" features also at small scales. For a recent and comprehensive overview of current issues and future prospects in understanding reactor neutrino spectra, see the contributions in [60].
In oscillation searches at reactor experiments, spectral uncertainties can be efficiently suppressed by comparing near (unoscillated) and far (oscillated) event spectra [61,62], as performed in [10][11][12]. In the context of JUNO, a concept for a high-resolution near detector was mentioned in [63] and further detailed in [64,65]. This concept has evolved into a full-fledged project, the Taishan Antineutrino Observatory (TAO) [66][67][68][69][70]. 1 TAO is expected to monitor the unoscillated spectrum emitted by one of the Taishan nuclear reactors, with a gain of about ×1/2 in energy resolution and ×30 in event statistics with respect to the oscillated spectrum at JUNO. Independently of neutrino oscillations, high-resolution spectral measurements at TAO will set unprecedented benchmarks [66] for research in nuclear fission physics [72][73][74] and for broader investigations of the neutrino-nuclear response in particle physics and astrophysics [75]. In general, progress in neutrino and nuclear physics, coupled with precision measurements at TAO, is expected to significantly constrain the range of neutrino spectral models to be used in future JUNO data analyses.
In this work we build upon our previous studies of precision oscillometry [76,77], but considering summation spectra with substructures and possible uncertainties. We use the publicly available toolkit Oklo [78,79] to generate ensambles of spectra within quoted errors on yields, branching ratios and endpoint energies for each decay. This toolkit, although currently not updated in terms of nuclear database inputs (taken as of 2015 [78]), is appropriate for our methodological purposes and numerical experiments. For simplicity, we shall assume that the underlying neutrino spectra are the same in TAO and JUNO. In reality, the former will closely monitor only one reactor core in Taishan, while the latter will detect a signal generated by several reactors in both Taishan and Yangjiang, with different fuel evolutions [66,71,80]. The related fuel corrections will require detailed information and modeling for each reactor, that are beyond the scope of this paper and will be studied elsewhere.
Our work can be divided in two main parts. In the first part (Secs. II and III) we discuss the formal relations between the TAO and JUNO energy spectra. We start by revisiting in detail the effects of resolution and recoil that, although well known in principle, are not always properly distinguished and implemented at the level of accuracy required by future measurements. Then we show that any observable energy spectrum of events in TAO can be mapped into a corresponding spectrum in JUNO by a proper convolution. In particular, we show that this mapping can be exactly performed in the hypothetical case of no oscillations (Sec. II) and can be very accurately generalized, via an ansatz, to the real case of neutrino oscillations (Sec. III). These results allow to predict the JUNO spectrum directly from a model for the observable event spectrum at the TAO detector, rather than from a model for the unobservable neutrino spectrum at the reactor source.
In the second part of the paper we perform quantitative studies of the mass-ordering sensitivity and precision oscillometry in JUNO, first by considering a single reference spectrum with substructures, and then by adding bundles of variant spectra to be constrained by TAO. In Sec. IV we revisit our previous analysis [77] including the reference Oklo spectrum, new priors for the oscillation parameters, and reduced error bands for smooth flux-shape and energyscale systematics. We confirm the accuracy of the mapping and discuss the impact of these new inputs. In Sec. V we generate bundles of Oklo variants around the previous reference spectrum. We perform a χ 2 analysis of variant spectra in JUNO, alternative to the Fourier analysis in [58], and highlight several statistical issues arising from sampling the nuclear input uncertainties. By varying all the known nuclear data inputs, we generate and analyze an ensemble of 10 5 spectra in JUNO, but find no reduction of the sensitivity to mass ordering, with or without TAO; we trace this unexpected result to subtle undersampling issues in the generated bundle of spectra. We repeat the analysis by constructing an equally numerous but "more densely sampled" bundle, and find a small reduction of the JUNO sensitivity, consistent with [58] and improved with the help of TAO. These results, based on "known" nuclear inputs, suggest some cautionary comments on parametrizations of "unknown" substructure uncertainties, as those considered in [57]. We also analyze the JUNO accuracy on the relevant mass-mixing parameters, which is found to be basically unaffected by fine-structure issues. Our results are summarized in Sec. VI.

II. MAPPING THE SPECTRUM FROM TAO TO JUNO WITHOUT OSCILLATIONS
Reactor neutrinos can be detected through the inverse beta decay (IBD) process ν e + p → e + + n followed by e + annihilation and delayed n capture. Using the notation of [76,77], we focus on two energy variables for IBD events: We also consider a less accurate approximation, dubbed as "mid-recoil," whereby the midpoint of the interval in Eq. (10) is taken as a proxy for E e [83], and the Jacobian is included, when passing from neutrino to positron energy spectra, to ensure event number conservation. A useful approximation for E mid e (and thus for J) is given in [83] as J(E) where ∆ E = m e + 0.783 MeV. This mid-recoil recipe captures well the average recoil shift but ignores its energy spread, which is definitely nonnegligible in TAO as shown below.

C. Resolution and recoil effects: comparison and combination
If both resolution and recoil effects were neglected, then Eqs. (6) and (9) would lead to the often-quoted approximation E vis = E − 0.783 MeV. Figure 1 (left panel) shows the recoil corrections to such relation as a function of neutrino energy E, in terms of deviations from unity [dashed line at 1 ≡ E vis /(E − 0.783 MeV)]. The gray area corresponds to the the one-sided energy deficit due to full recoil effects [Eq. (10)], while the solid line marks the mid-recoil approximation [Eq. (12)]. Notice that, at high reactor neutrino energies, the visible event energy is both shifted and smeared out at the percent level. In Fig. 1 (right panel) we show the fractional energy spread ∆E vis /E vis due to recoil and resolution, separately. In particular, ∆E vis is shown as ±(E 2 − E 1 )/2 for recoil (dark gray), as ±σ T for TAO (gray band) and as ±σ J for JUNO (light gray). Recoil and resolution effects in TAO appear to be of comparable size, and none of them can be neglected in accurate spectral analyses, especially in view of their impact on the observability of substructures. As shown in [76], the combination of the resolution and recoil effects is fully encoded in an energy resolution function R X that connects the relevant energies E vis and E, as obtained by convolving the gaussian distribution r X in Eq. (7) with the top-hat distribution t in Eq. (11), where σ X = σ X (E vis ), while the dependence on E comes from E 1,2 = E 1,2 (E) that we take from the full calculation in [82]; see [76] for further details, including the adopted convention for the error function (erf). The observable energy spectra S X of IBD events in TAO and JUNO (in the absence of oscillations) can then be computed by convolving the neutrino spectrum S ν in Eq. (3) with the above resolution function, where E T = 1.806 MeV is the IBD ν energy threshold and N X is a normalization factor. Figure 2 shows in the left panel the neutrino energy spectrum S ν , as obtained with default nuclear input parameters for the ν flux Φ from the Oklo toolkit [78], times the cross section from [82]. The TAO visible energy spectrum S T is shown in the right panel, including recoil and energy resolution effects. Spectra are normalized to the same area (in arbitrary units) to facilitate comparison in shape. It can be seen that spectral substructures in S ν (sawtooth and step-like features) are smeared out in S T , but still partly visible. Such substructures would no longer be visible in JUNO (not shown). Figure 3 shows the ratio S T /S J of the unoscillated energy spectra in TAO and JUNO (arbitrarily normalized to the same area). From left to right, the numerator S T is calculated with progressive inclusion of recoil and resolution effects, while the denominator S J always includes all such effects. In particular, the first three panels assume perfect energy resolution in TAO (σ T = 0), with increasingly accurate treatments for nucleon recoil. In the first plot (no recoil) the spectral ratio shows evident substructures and a high-energy excess (spectral tilt), due to neglected energy recoil losses that also bias the substructure peak positions by up to 1% (not visible by eye). In the second plot (mid-recoil approximation including the Jacobian) the average energy losses are accounted for, the shift disappears, and the subtructures are correctly aligned in energy. In the third plot (full recoil treatment) the inclusion of the recoil  The denominator S J always include full recoil and resolution effects. The numerator S T includes recoil and resolution effects in progression from left to right. In the rightmost plot, the spectral ratio substructures are compared with the ±1 statistical error band in TAO, assuming 3 × 10 6 events and 40 keV bin width. energy spread suppresses the finest spectral structures and, at high energy, reduces their amplitudes by a factor of ∼ 2. Finally, further suppression of fine structure features (and another amplitude reduction by a factor of ∼ 2 or more) is due to the inclusion of the finite TAO resolution width σ T from Eq. (8) in the rightmost panel. In this panel we also show the ±1σ error band in TAO assuming 3 × 10 6 events, i.e., ∼ 30 times the statistics expected in JUNO [66] in the presence of oscillations for about 5 years (that amounts to ∼ 100, 000 events [76]). The statistical band depends on the bin width, here taken as 40 keV (25 bins per MeV interval) in order to cover the most prominent substructures within a few bins at least. It can be seen that a handful of fine-structure features reaches the ∼ 1σ level in amplitude, allowing TAO to probe spectral models with different predicted substructures (see also [66]). We shall discuss some statistical issues concerning the model selectivity of TAO in Sec. V.
Summarizing, resolution and recoil effects in the TAO energy spectrum are of comparable size and should be carefully implemented, in order to avoid energy biases and unrealistic amplitudes for fine-structure spectral features. Resolution effects produce a gaussian smearing (whose width decreases with increasing energy), while recoil effects produce an energy shift plus a top-hat smearing (whose width increases with increasing energy). Their combination (convolution) leads to an analytical expression for the energy resolution function [76] as in Eq. (16), that can be usefully applied to the calculation of both TAO and JUNO spectra.
Finally we mention that, in principle, the impact of recoil effects may be reduced by directional information in the final state of IBD events, see [84] for a recent proposal in the context of TAO. We do not explore this option hereafter, but surmise that constraining recoil effects amounts to replace the function t in Eq. (11) with another one (t ) having smaller variance, possibly leading to an analytical result as in Eq. (16) if the parameterization of t is simple.

D. Mapping the spectrum from near to far
The resolution function R J in Eq. (16) for JUNO is obtained by convolving a gaussian r J having a variance σ 2 J with a top-hat function t. In turn, r J can be thought as the convolution of two gaussians r T and r D with variances given, respectively, by σ 2 T (as in TAO) and by that is, the difference between the energy resolution variances in JUNO and TAO. Then, through convolutions, one gets an exact mapping from TAO to JUNO (unoscillated) spectra as follows: where normalization factors N X have been dropped for simplicity, and r has the same functional form as in Eq. (7), with E e + m e replaced by E vis .
This analytical result has a simple physical interpretation: The JUNO unoscillated spectrum in visible energy (S J ) can be obtained from the TAO spectrum (S T ) by applying an extra gaussian smearing with variance σ 2 D , equal to the difference of variances in JUNO (σ 2 J ) and TAO (σ 2 T ). In doing so, recoil effects remain correctly implemented in both TAO and JUNO.
Note that Eq. (19) directly relates the observable event spectra S T and S J , without using the unobservable neutrino spectrum S ν . This represents an advantage in terms of nuclear physics modeling: Constructing a model for S T (compatible with future TAO data) will generally be less demanding than building a complete model for S ν , since the former will exhibit only a few surviving substructures to be properly described via summation.
A final comment is in order. As stated in Sec. I, we are assuming the the TAO and JUNO spectra are generated by the same underlying ν spectrum S ν . However, JUNO will collect a neutrino flux also from reactor cores different from the one monitored by TAO, leading to fuel-component differences in the reference S ν and to corrections to the ideal case in Eq. (19). Fuel evolution issues and related spectral effects in TAO versus JUNO are beyond the scope of this investigation, and will be treated in a future work; see [66,71,80] for useful considerations in this context.

III. MAPPING THE SPECTRUM FROM TAO TO JUNO WITH OSCILLATIONS
In this Section we generalize the TAO → JUNO mapping of Eq. (19) in the presence of oscillations, characterized by a ν e survival probability P ee (E). An obstacle to this goal is that the integrand S ν (E) gets replaced by the product S ν · P ee in JUNO, and that the convolution of a product is not the product of convolutions, as also noted in [80].
However, after reviewing the functional form of P ee , we propose an ansatz that, to a very good approximation, overcomes this problem. We shall generalize Eq. (19) by including an effective probability P eff ee , expressed in terms of observable spectra S X and visible energy E vis , that bypasses any prior knowledge of the (unobservable) neutrino energy spectrum S ν (E). We shall then discuss the validity of this ansatz, and use it in an updated analysis of the JUNO sensitivity to the neutrino mass ordering and to precision oscillometry.

A. Oscillation probability in terms of neutrino energy
In this subsection we describe the survival probability P ee (E), largely following [76] to which we refer the reader for details and references. In general, P ee (E) in JUNO depends on several parameters, where δm 2 = m 2 2 − m 2 1 and ∆m 2 = |m 2 3 − (m 2 1 + m 2 2 )/2| > 0 are the squared mass splitting parameters, α = ±1 distinguishes the mass ordering (normal or inverted), θ 12 and θ 13 are the mixing angles, N e is the electron density in matter, and {w n , L n } characterizes the set of reactors, each contributing to the total flux with fractional weight w n ( n w n = 1) at distance L n , under the assumption of identical fuel components.
Useful derived parameters are where c 12 = cos θ 12 and s 12 = sin θ 12 , and where L = n w n L n is the average baseline. Matter effects in JUNO depend on the ratio µ = (2 √ 2G F N e E)/δm 2 and lead to an effective mass-mixing pair (δ,θ 12 ) given byδ δ(1 + µ cos 2θ 12 ) and sin 2θ 12 sin 2θ 12 (1 − µ cos 2θ 12 ) at first order in the small parameter µ [76] (see also [86,87]). The dependence of P ee on its parameters can then be simply expressed as encodes (δ,θ 12 ) matter effects, while w 1−2∆ 2 ee n w n (1−L n /L) 2 is a damping factor due to the spread of baselines L n , and ϕ is the interference phase directly related to mass ordering [85]. An accurate empirical parameterization for ϕ is given by [ where P reads as in Eq. (24) but with vacuum mass-mixing values (δ, θ 12 ). 3 For the oscillation parameters in P ee [Eq. (23)] we assume the following priors (central values and ±1σ, after symmetrizing errors and averaging NO-IO differences) from the global analysis in [88]: Determining the mass ordering in JUNO amounts to prove that in P ee , besides the oscillation phase 2∆ ee , there is and extra an interference phase ϕ endowed with a definite sign (α = ±1) and not scaling as 1/E, otherwise it would be absorbed into a shift of ∆m 2 ee [89]; equivalently, one should find evidence for a non-constant ratio ϕ/2∆ ee . It has been pointed out [90] that energy calibration errors at (sub)percent level may (partly) mimic ϕ/2∆ ee = const [15,76,77]; in this context, future evidence for some substructures emerging in TAO spectrum, located at the energies predicted by nuclear summation models, may help the overall calibration of the reference spectrum to be projected from TAO to JUNO (provided that JUNO is also accurately calibrated in energy). Correct implementation of recoil effects, in both TAO and JUNO, remains mandatory to avoid energy biases at comparable (sub)percent levels. Figure 4 shows the function P ee (left panel) and the ratio ϕ/2∆ ee (right panel) as a function of energy. The solid curves and gray bands correspond, respectively, to the central values and to the envelopes of ≤ 1σ variations for the oscillation parameters. Normal ordering is assumed. The smallness of ϕ/2∆ ee illustrates the challenges of mass ordering determination at MBL reactors. Note the relatively high values of ϕ/2∆ ee for E ∼ 3 MeV may fractionally change by up to ±8% within the gray band, and that for similar energies P ee (and thus the IBD event rate) may also change by up to ±12%. Therefore, variations of the oscillation parameters within their current global-fit errors can appreciably affect the prospective mass ordering sensitivity in JUNO, as discussed later.
B. Ansatz: effective probability in terms of visible energy Given the probability P ee (E), the oscillated spectrum at JUNO (including resolution and recoil effects, and up to a normalization factor) is Our goal is to obtain such S J by mapping the TAO spectrum S T , in a form analogous to Eq. (19), where P eff ee should act as an effective oscillation probability, expressed in terms of the measured visible energy rather than the unobservable neutrino energy. This problem is exactly solved by imposing, in the kernel of Eq. (31), that namely, by defining P eff ee as follows (with a change E vis → E vis in the dummy variable): which represents the weighted average of P ee over the neutrino spectrum (S ν ) times the TAO energy resolution function (R T ). In a sense, P eff ee is a smeared version of P ee , averaged over S ν variations on an energy scale set by σ T . However, this formally exact solution is not satisfactory, as it requires the knowledge of the unobservable neutrino spectrum S ν . We make then the following ansatz, that replaces the unobservable S ν with its closest observable proxy, namely S T : Within the integral kernels of Eq.
The effect of the Jacobian in the above formula is rather small numerically, since J(E) changes slowly with E (if it were constant, it would be canceled in the ratio); we keep it for the sake of completeness. Summarizing, our ansatz for the mapping S T → S J (including oscillations) consists in calculating an effective JUNO spectrum S eff J (E vis ) from the observable TAO spectrum S T (E vis ) as via the effective probability In the limit of no oscillations (P ee = 1 = P eff ee ), Eq. (35) reproduces the exact result in Eq. (30). that this recipe can approximately capture the local smearing of P ee implicit in Eq. (33) without introducing energy biases, as the average recoil effects are accounted for by the mid-recoil approximation. Of course, the replacement of S ν (E) with the proxy S T E mid e (E) + m e introduces an extra smearing associated to the latter spectrum, which is absent in the former. This artifact may be expected to have marginal effects in the final S J , since the smearing in JUNO acts on an energy scale σ J > σ T . Ultimately, the validity of our ansatz relies on numerical tests. Figure 5 shows the ratio of the JUNO spectra calculated with the ansatz [S eff J from Eqs. (35) and (36)] and without the ansatz [S J from Eq. (30)]. The underlying neutrino spectrum S ν is taken as the reference Oklo spectrum in Fig. 1  (left panel). The solid line and gray band refer, respectively, to central values of the oscillation parameters and to their ±3σ variations (applied to both S eff J and S J at the same time). The ansatz provides numerically accurate results at the level of few ×10 −4 , except in the high-energy tail where it reaches a permill level, that is anyway insignificant as compared with other sources of uncertainties (both statistical and systematic) in JUNO, as also discussed below. Finally, we have tested that the same excellent accuracy in Fig. 5 is reached by replacing the reference spectrum with variant spectra, as obtained from the Oklo toolkit by changing the nuclear inputs within their uncertainties (not shown).

IV. NEUTRINO OSCILLOMETRY IN JUNO: SINGLE SPECTRUM
We present and discuss a prospective analysis of JUNO in terms of sensitivity to mass ordering and of precision determination of oscillation parameters, building upon our previous work [77]. Here we use a single input spectrum, namely, the reference Oklo one as shown in the previous Sections. Bundles of variant spectra and their effects will be considered in the next Section. The main purpose of this updated analysis is to further test the previous ansatz and to discuss the impact of changes in the reference oscillation parameters and other systematics. TAO does not play a specific role herein, except for providing a reference spectrum S T for the S T → S J mapping, when the ansatz is used.  A. Ingredients of the analysis Figure 6 shows the observable JUNO spectrum S J expected in the presence of oscillations from the Taishan and Yangjiang reactor sources (dashed red line) plus the background components expected from farther reactors (blue dotted line) and U+Th geoneutrinos (green solid line). 4 The total spectrum (black solid line) is endowed with a gray band, representing the envelope of variations of the oscillation parameters within their prior 1σ ranges. All curves refer to 5 years of data taking (∼ 10 5 JUNO events), assuming the same normalization factors for the various components as discussed in [77], to which we refer the reader for details not repeated herein.
We focus here on the inputs that differ from [77]. The central values (and to some extent the errors) of the oscillation parameters in Eqs. (26)-(29) have changed, in particular for the mass splittings (about +1σ for ∆m 2 ee and −1σ for δm 2 ee ). Concerning Φ, we use the reference neutrino flux from the Oklo toolkit, corresponding to the neutrino spectrum S ν = Φσ ν in Fig. 1 (left panel). Note that, in this Section, we do not attach uncertainties to the fine structures of Φ, that will be separately addressed in Sec. V. However, we do include large-scale (smooth) uncertainties of the flux shape in the form of polynomial deviations Φ /Φ, as well as energy-scale systematics in the form of polynomial deviations E /E, adopting the same methodology as in [77] but with narrower error bands. Figure 7 shows our default ±1σ error bands for Φ /Φ and E /E variations (left and right panels, respectively). Here we reduce the width of the Φ /Φ band to 2/3 of the previously adopted one in [77], because: (1) at low energy, the normalization error (that sets the lower limit to the width) has been reduced from ∼ 2.3% [77] to ∼ 1.5% [28]; at high energies, prospective analyses of the flux shape reconstruction in TAO [66] give reasons for moderate optimism. The E /E error band is taken from [91] (see Fig. 18 therein), with an appreciable reduction (roughly by a factor 1/2) with respect to [77].

B. Sensitivity to mass ordering
Following [77], we perform a least-squares analysis of the JUNO sensitivity to mass ordering, up to 10 years of data taking. We remind that our complete χ 2 function for JUNO is defined as where: the first term includes statistical errors only; the second term includes penalties for variations of the oscillation parameters, governed by the priors in Rqs. (26)-(29); the third term contains normalization errors for the geo-ν Th and U fluxes, normalization and (polynomial) shape systematics for the reactor fluxes, and (polynomial) energy-scale systematics. The second and third term contain up to N sys = 18 systematics, treated as nuisance parameters that are marginalized away in the χ 2 J minimization [77]. The analysis is performed by progressively including such nuisance parameters: (1) oscillation parameters and normalizations (osc. + norm.), N sys = 7; (2) plus energy scale variations, N sys = 13; (3) plus flux shape variations, N sys = 18. Normal (inverted) ordering is assumed as true (test) hypothesis.  Note that systematics do not seem to saturate the sensitivity to mass ordering even with 10-years data. Also note that this sensitivity is reduced more by energy-scale uncertainties than by flux-shape ones. Therefore, it will be important to ensure that the energy calibration in JUNO can achieve the same (or better) level of accuracy reached in [28].
It is useful to compare this Fig. 8 with the analogous Fig. 7 in [77]. It turns out that the curve N σ (T ) including the full set of systematics (red curve) is almost unaltered, despite the previously discussed reduction in energy scale and flux shape systematics. The (surprising to us) explanation is that the benefits of this error reduction happen to be accidentally compensated by "unlucky" changes in the central values of the oscillation parameters. Further understanding can be gained by focusing on the case with "osc. + norm." errors only (black curve in Fig. 8) for a fixed live time T = 5 years. plane (δm 2 , sin 2 θ 12 ) for fixed (∆m 2 ee , sin 2 θ 13 ), and viceversa in the right panel. The coordinates span the ±2σ ranges for the mass splitting and ±1σ ranges for the mixing angles, in the units of Eqs. (26)- (29). The ∆χ 2 value increases noticeably by increasing δm 2 or by decreasing ∆m 2 ; in other words, the mass ordering test in JUNO improves when the ratio ρ = δm 2 /∆m 2 ee increases, even if by small amounts (conversely, the mass ordering would become eventually untestable for vanishing ρ). Note that a (more modest) increase of ∆χ 2 is also obtained by increasing either sin 2 θ 12 or sin 2 θ 13 and thus the oscillation amplitude(s), as it can be generally expected in oscillation searches.
It turns out that, with respect to [77], the central values of all four oscillation parameters in Eqs. (26)-(29) have accidentally changed in "unlucky" directions, lowering ∆χ 2 by about 3.5 units for the case of "osc. + norm." uncertainties. As anticipated, for the analysis including all the uncertainties, this drop is almost exactly compensated (once more, accidentally) by the reduction of energy-scale and flux-shape systematics. Similar results have been obtained in the case where the true ordering is inverted and NO is tested (not shown). In conclusion, the JUNO rejection of the wrong mass ordering depends, in a nonnegligible way, on the central values of the oscillation parameters.

C. Accuracy of oscillation parameters
Eventually, at least three oscillation parameters (δm 2 , ∆m 2 ee , θ 12 ) will be very precisely measured by JUNO itself. Figure 10 shows the time evolution (in JUNO) of the fractional accuracy σ p /p for each of six parameters p, namely, from top to bottom: the two mass splittings, the two mixing angles, and the U and Th geoneutrino flux normalizations (f U and f Th ). For each parameter, it is understood that the others are marginalized away in the analysis. At T = 0, color code is the same as in Fig. 8, and the abscissa also scales as √ T . From top to bottom, the results refer to the squared mass splittings δm 2 and ∆m 2 ee , the mixing angles sin 2 θ 12 and sin 2 θ 13 , and the geoneutrino flux normalization factors f U and f Th . The same results are obtained by using the S T → S J mapping ansatz (not shown). the oscillation parameter errors are set by Eqs. (26)- (29), while for the geo-ν fluxes we assume the priors in [77], f U = 1±0.20 and f Th = 1±0.27. The color sequence for the curves (red, blue and black for growing sets of systematics) is the same as in Fig. 8. After a live time of 5 years, the accuracy of (δm 2 , ∆m 2 ee , θ 12 ) will improve by factors of about (6, 6, 4), respectively-or better, if some systematics che be further reduced. For the pair (δm 2 , ∆m 2 ee ), that governs the "slow" oscillations in the JUNO spectrum, flux-shape uncertainties are more important than energy-scale ones, and viceversa for ∆m 2 ee that governs the "fast" oscillations. A moderate reduction of the prior errors will be obtained for geoneutrino fluxes, with little dependence on systematics. Concerning θ 13 , the current experimental error will only be marginally improved. Finally, we have repeated the analysis by using the S T → S J mapping ansatz, obtaining the same results with insignificant deviations (not shown).

V. NEUTRINO OSCILLOMETRY IN JUNO: ENSEMBLES OF SPECTRA
Summation calculations of reactor spectra have come a long way since the pioneering works [92][93][94][95]. Modern realizations are based on thousands of nuclear input data on decay yields Y i , endpoints Q j and branching ratios b k , together with their uncertainties and possible covariances [56,96]. However, as mentioned in the Introduction, even the most refined summation spectra do not match well current reactor data, suggesting that some nuclear (experimental or theoretical) ingredients may be missing. Significant work is still needed to reach consensus on satisfactory spectra with realistic uncertainties and correlations [60,97], with TAO providing important benchmarks in the future [71].
With all these caveats, we perform an exploratory analysis of the effect of "known" nuclear input uncertainties on the spectral substructures through the Oklo toolkit [78]. We remind that the Oklo code contains (4306, 6609, 6804) values for (Y i , Q j , b k ), respectively, for a total of N d = 17, 719 input data, together with their quoted uncertainties (taken as uncorrelated). These huge numbers prevent usual χ 2 analyses of variant spectra, in terms of marginalization over nuisance parameters. Alternatively, we generate ensembles of N neutrino spectra {S n ν (E)} n=1,...,N , by randomly varying all or some nuclear inputs within their uncertainties. We also compute the associated TAO spectra {S n T (E vis )}, that are then mapped to obtain JUNO spectra {S n J (E vis )} (where we drop the superscript "eff" for simplicity). We test how these variants affect the JUNO oscillation analysis, and how well they can be distinguished by TAO, by scanning appropriate χ 2 functions over the whole spectral set(s). We recover, through an independent χ 2 analysis, the results obtained in [58] through a Fourier analysis, namely, that "known" substructure uncertainties do not appear to pose a threat to precision oscillometry in JUNO. However, the quantification of this result is not trivial, and some subtle problems in the statistical analysis will be highlighted. We shall also comment on the issue of possible "unknown" small-scale uncertainties, as raised e.g. in [57] and in [59,80].
A. Changing all nuclear input uncertainties: spectrum metric and (under)sampling issues In our first exercise with spectral variants, we have generated an ensemble of N = 10 5 neutrino spectra S n ν (and associated TAO spectra S n T ) by N extractions of random values s n i for all the i = 1, . . . , N d inputs at the same time, assuming uncorrelated gaussian distributions for the quoted uncertainties σ i . At each extraction, branching ratios for each decay are renormalized by an overall factor to ensure unitarity ( k b k = 1). All variant spectra S n T are normalized to the same area as the reference spectrum S T , in order to emphasize shape variations. Figure 11 shows S T (solid line) with its statistical errors (dark gray band), assuming 3 × 10 6 IBD events in TAO, and 40 keV bins. Also shown is the envelope of all the S n T variant spectra (light gray band), and a few individual variants (very light gray curves). All spectra are divided by the unoscillated JUNO spectrum S J , analogously to Fig. 3. Since the light gray band is rather large, one may expect that at least some spectral variants within the envelope can play a role in the TAO and JUNO data analyses. The surprising outcome is that only the reference S T matters in our exercise, for subtle reasons that we could not anticipate.
A first issue is how to define a χ 2 S metric within the {S n T } envelope, so that 68% (95%) of the spectra fall within a properly defined 1σ (2σ) band etc. (with N σ = χ 2 S ) around the reference spectrum S T . Note that each spectrum S n T is endowed with a χ 2 n value that measures its statistical distance from S T (having χ 2 = 0 by definition) in terms of nuclear input uncertainties. For N d 1, the distribution of χ 2 n values (not shown) can be approximated by a gaussian centered at N d and with variance 2N d [98], effectively starting at 0 (corresponding to S T ) rather than −∞. Since S T sits in the tail rather than at the peak, this distribution does not directly provide a good metric. However one can construct a proper metric χ 2 S by integrating this χ 2 n distribution from zero up to the fractional area corresponding to the desired N σ level. As a result (proof omitted), each spectrum S n T is endowed with a new χ 2 S,n value given by: where erf −1 is the inverse error function, and χ 2 S = 0 is recovered for S T in the limit N d 1. It turns out that if the bands corresponding, e.g., to χ 2 S,n ≤ 1, 2, and 3 were plotted in Fig. 11, they would be insignificantly smaller than the light gray envelope of all the spectra. In other words, by taking spectra with increasingly high χ 2 S,n (or equivalently χ 2 n ), more variant spectral shapes become possible within the band, while the typical subtructure amplitudes remain constant and their envelope is not enlarged.
These results suggest caution in parametrizing variant spectra as in [57], namely, by breaking down the envelope in bins and computing uncorrelated standard deviations in each bin, for two reasons: (1) the amplitude of deviations does not scale with N σ ; (2) by binning, the detailed information about which shapes are (not) allowed by nuclear uncertainties is completely lost; in particular, the loss of point-to-point correlations permits more shapes than would be allowed by nuclear inputs only. In doing so, "known" uncertainties are partly replaced by "unknown" ones, allowing substructure amplitudes and shapes beyond those pertaining to compiled nuclear inputs.
A second issue concerns the fraction of spectra {S n T } that survives the comparison with prospective TAO data. We consider a simplified χ 2 analysis for TAO, where each spectrum S n T is compared with the reference one S T in terms of statistical errors, plus one nuisance normalization parameter λ (S n T → λS n T , assuming σ λ /λ = 1.5 × 10 −2 ), in addition to χ 2 S,n that embeds nuclear errors: where for χ 2 stat,n we adopt the limit of infinite bins [76,77,99], that provides a very good approximation to the binned case. Within the ensemble {S n T }, the fraction of spectra allowed at N σ by TAO data (defined by χ 2 TAO,n ≤ N 2 σ ) is a function of the TAO exposure. With ∼ 3 × 10 6 IBD events expected in TAO after ∼ 5 years, we unexpectedly find that none of the 10 5 synthetic spectra survives, even at N σ = 3 level: they are all rejected with respect to the reference spectrum S T . It turns out that the good TAO energy resolution is sufficient to distinguish spectra S n T that differ from S T by a few substructures, even with much less than ∼ 3 × 10 6 events. Figure 12 shows how well TAO selects spectra in the ensemble {S n T }, as a function of the total number of collected IBD events. The three curves represent the fraction of spectra that survives at N σ = 1, 2 and 3. For no TAO events these fraction correspond by construction, respectively, to 0.68, 0.95 and 0.997. By increasing the number of TAO events, these fractions drop rapidly. When the surviving fractions drop below 10 −4 (not shown), the curves break down because only a handful of the 10 5 spectra-and ultimately none of them-is allowed; this happens well below 3 × 10 6 events in TAO. The results in Fig. 12 suggest that, when all the nuclear input uncertainties (N d = 17, 719) are randomly varied, generating 10 5 synthetic spectra by random extractions is not enough to densely sample the ∞ N d -dimensional set of possible variant spectra: orders of magnitude more extractions would be needed to obtain a few spectra S n T really close to S T within statistical uncertainties. A third statistical issue, connected with the last one just discussed, concerns the JUNO sensitivity to mass ordering. We have repeated the prospective JUNO data analysis in Sec. IV by mapping the spectral ensemble {S n T } → {S n J } for any set of oscillation parameters. In particular, assuming NO and the reference S J as the true hypothesis, we have tested the wrong IO not only via S J but also scanning the 10 5 spectra S n J (with and without adding the term χ 2 S,n ). 5 We have found no reduction of the sensitivity to the mass ordering as compared with Fig. 8. 6 These results qualitatively agree with those in [58] (where a Fourier spectral analysis found that substructures played a little role) but are unexpectedly stronger: none of the test spectra induces any sensitivity reduction in JUNO. In addition, we find that also the precision determination of several parameters p as in Fig. 10 remains unaltered. Once more, we surmise that the ensemble of 10 5 spectra is not dense enough to sample shape variations very close to the reference one. In order to overcome these issues we construct and test a denser ensemble below.

B. Changing only some nuclear input uncertainties: Suggestions for possible parametrizations
We have constructed an alternative ensemble of 10 5 spectra {S n T } with substructures closer to the reference spectrum S T as follows: at each of 10 5 extractions, we have randomly chosen a subset of only N d = 10 2 nuclear input uncertainties (out of N d = 17, 719) to be varied. Figure 13 is analogous to Fig. 11 but shows the envelope of such new spectra, which is narrower and closer to S T by construction. Also in this case, by ranking variant spectra with a χ 2 S,n as in Eq. (39) (with N d replaced by N d ), the N σ bands would be only marginally smaller than the envelope, confirming that substructure amplitudes do not scale with N σ .  Figure 14 is analogous to Fig. 12 but with the new ensemble of spectra. In this case, O(10 6 ) TAO events are required to start reducing the fractions of spectra allowed at N σ . For 3 × 10 6 events, using Eq. (40), we find that the envelope of spectra surviving at 1σ is as shown in Fig. 15. The envelopes at 2σ and 3σ (not shown) are about a factor of 2 and 3 larger than the light gray band in Fig. 15, suggesting that the fit to TAO data tends to linearize the scaling of allowed substructure amplitudes with N σ .
Concerning the JUNO sensitivity to mass ordering, we now find a slight reduction of ∆χ 2 (IO − NO), amounting to −0.4 when scanning over the whole new set of {S n J }; this reduction is halved to −0.2 when this set is reduced by TAO via Eq. (40). These relatively small effects, derived through a χ 2 analysis, agree with the Fourier-analysis findings of [58]: variant spectral substructures appear to play a little role in the JUNO sensitivity to mass ordering, as     FIG. 14: Spectral ensembles in TAO. As in Fig. 12, but considering only a random subset of 100 nuclear uncertainties.  Fig. 13, but with the light gray band representing the envelope of spectra allowed at 1σ by TAO after collecting 3 × 10 6 events.
far as known nuclear uncertainties are concerned. 7 The role is even more marginal with the help of TAO. Of course, if all substructures shapes were hypothetically allowed, including oscillatory ones appropriately tuned to "undo" the IO-NO probability differences, then the sensitivity reduction would be higher [57], at the price of introducing ad hoc "unknown" errors, not belonging to those parametrized in nuclear databases.
From the exercises in this Subsection and in the previous one we learn that, once TAO spectral data and an associated reference model spectrum S T (E vis ) will be available, it will be important to sample very densely the functional neighborhood of such spectrum, in order to study the residual effects of allowed variant spectra in JUNO. Brute force variations of all the O(10 5 ) nuclear decay parameters may lead to undersampling issues in this context. Reduction to a limited number of nuclear error sources appears to be a better strategy. While we have arbitrarily limited this number to 10 2 random error sources, future studies may motivate on more physical grounds a limited subset of nuclear errors (plus possible covariances), related to the decays producing the most pronounced substructures in TAO. If the nuclear physics of reactor neutrino spectra will not be well understood even in in the TAO era, these "known" error sources may be cautiously supplemented (but not replaced) by some extra errors for "unknown" substructures.

VI. SUMMARY AND CONCLUSIONS
The next-generation, medium baseline reactor neutrino experiment JUNO (in construction) is planned to probe the full pattern of ν e disappearance for L/E ∼ (few MeV)/(53 km), including the precision measurements of oscillations induced by the (δm 2 , θ 1 2) and (∆m 2 , θ 13 ) mass-mixing pairs, and their interference effects governed by the neutrino mass ordering, namely sign(δm 2 /∆m 2 ). The supplementary detector TAO is expected to monitor the unoscillated flux close to one reactor core, with about a factor ×2 improvement in energy resolution and with ×30 more events than in JUNO.
In this work we have studied the relations between the observable event spectra in TAO (S T ) and JUNO (S J ), in the simplifying assumption that they are generated by the same unobservable neutrino spectrum (S ν ), including fine-structure features as emerging in summation calculations. We have used the publicly available Oklo toolkit [78] to generate a reference spectrum S ν , as well as a number of variants S n ν corresponding to changes in the (thousands) of nuclear inputs describing fission yields, branching ratios and endpoint energies. Our methodology can be applied to more updated nuclear databases, which are currently being developed and endowed with preliminary error covariances (not included in this work). After reviewing in detail the different and non-negligible effects of energy resolution and nucleon recoil on the observable spectra, we have shown that a model spectrum S T at TAO site can be mapped into a corresponding spectrum S J at JUNO via well-defined convolutions, without using the (more detailed) information contained in the source neutrino spectrum S ν . The mapping S T → S J is exact in the hypothetical case of no oscillations, and can be generalized with excellent accuracy to the real case with oscillations, via an ansatz on the effective disappearance probability. The prospective χ 2 analysis of JUNO data confirms the validity of the mapping, and allows to discuss the impact of uncertainties related to oscillation parameters, energy-scale and flux-shape systematics.
We have also analyzed the effect of known nuclear input uncertainties, by generating bundles of variant spectra with the Oklo toolkit. We highlight several statistical issues arising from sampling a large number of variable inputs. We find that the bundles must densely sample the neighbourhood of the reference spectrum, in order to produce a detectable effect on the JUNO χ 2 function in numerical experiments. In this case (realized by sampling only a random subset of nuclear uncertainties), the effect turns out to be small (in agreement with [58]), and can be further reduced by adding TAO constraints. These results, based on "known" nuclear inputs, also suggest some cautionary comments on parametrizations of "unknown" substructure uncertainties, in terms of variances of binned bundles.
We have argued that, when TAO data will be available, an optimal strategy to deal with small-scale spectral shape uncertainties will be to focus on a few prominent visible substructures and related nuisance parameters, in order to build a dense ensemble of TAO spectral variants, to be mapped in JUNO and compared with its data. Optimal constructions for such variant ensembles, possibly with covariances of known errors and with allowance for extra unknown errors, as well as for corrections due to different fuel components in the TAO and JUNO sources, are left to future studies. We conclude by observing that, after only two decades from the discovery of neutrino oscillations, the JUNO and TAO projects are bringing the discussion of precision oscillometry to an unprecedented level of details, whose deeper understanding will require further advances at the junction of neutrino and nuclear physics.