Counting Reactor Antineutrinos: Huber, Mueller, and the Daya Bay Audit

A modern reactor antineutrino experiment compares two numbers: the number of inverse-beta-decay events its detector actually records, and the number it expects to record based on the predicted flux from the reactor cores it observes. The ratio is what carries the physics. If the experiment is at a short baseline where oscillation has not had time to set in, the ratio should be one; deviations point to flux miscalculation, sterile oscillation, or both. If the experiment is at an oscillation-baseline like Daya Bay’s 1.5-km far detectors, the ratio is set by the ${\theta }_{13}$ oscillation parameter and the predicted unoscillated flux serves as the reference.

The whole framework rests on knowing the predicted unoscillated flux correctly. This is harder than it sounds. A power reactor is a sprawling source: hundreds of fission fragments contributing simultaneously, fuel composition shifting over months as ${}^{235}$U depletes and plutonium isotopes accumulate, and antineutrino spectra that have to be inferred from beta-spectra measurements made in the 1980s with technology that doesn’t quite match what modern experiments need. The standard reference — the Huber-Mueller model — has held up well in its overall shape but produces a curious 6% normalisation deficit relative to the world’s reactor measurements (the “reactor antineutrino anomaly”) and a bump in the 5-MeV region that no model predicts. Daya Bay’s 2019 isotope-decomposition measurement provided the cleanest test of where the discrepancy lives.

This post is about how reactor antineutrino predictions are actually built, what the dominant uncertainties are, and how Daya Bay’s per-isotope measurement helped pinpoint where the models go wrong.

The basic decomposition

The predicted antineutrino flux at any detector position is

$\Phi (E_{\nu })={\sum }_i{N}_{f,i}{S}_i(E_{\nu })$

where $N_{f,i}$ is the rate of fissions of isotope $i$ in the reactor core and $S_i(E_{\nu })$ is the antineutrino energy spectrum per fission of isotope $i$, summed over all fission-fragment beta decays. The four major fissioning isotopes in a typical light-water reactor are:

• ${}^{235}$U (the principal fissile fuel, depleted as the cycle progresses)
• ${}^{238}$U (fast-fission and breeder, roughly constant)
• ${}^{239}$Pu (bred from ${}^{238}$U over the cycle, becomes substantial late in the cycle)
• ${}^{241}$Pu (bred at lower rates, contributes a few percent late in the cycle)

The fission rates $N_{f,i}$ are calculated by reactor physicists using core-simulation codes that model the neutron transport and the burnup over time. These calculations are quite reliable — typical uncertainty in the fission fraction is at the 1-2% level for each isotope. The neutrino-physics contribution to the prediction is the spectrum $S_i(E_{\nu })$, which is where the modelling difficulty lies.

The conversion approach

The original and historically dominant approach to building $S_i$ is the conversion approach: take a measured cumulative beta spectrum from a sample of the fissioning isotope, then invert it into the corresponding antineutrino spectrum using kinematic relations from beta decay.

The cumulative beta spectra were measured in the 1980s at the Institut Laue-Langevin (ILL) high-flux reactor in Grenoble. Thin foils of pure ${}^{235}$U, ${}^{239}$Pu, and ${}^{241}$Pu were irradiated with thermal neutrons in the ILL core, and the beta spectra emitted by the resulting fission-fragment mixture were measured with a magnetic spectrometer. These measurements provide the total beta spectrum from all the unstable fragments after their respective decay chains, exactly the same way the antineutrinos come out of an operating reactor — just observed through their charged-lepton siblings.

Converting the measured beta spectrum to the antineutrino spectrum exploits a simple kinematic relation: in a beta decay with Q-value $Q$, the beta and the antineutrino share the available energy, so the beta spectrum at energy $E_{\beta }$ is in one-to-one correspondence with the antineutrino spectrum at energy $E_{\nu }=Q-E_{\beta }$. The problem is that the ILL cumulative spectrum sums over hundreds of fragments with different Q-values, and the conversion has to handle the convolution carefully. The standard approach fits the measured cumulative spectrum to a sum of about 30 “virtual” branches, each with a distinct Q-value and weighting, and then maps each branch into the corresponding antineutrino contribution.

Patrick Huber’s 2011 reanalysis applied modern nuclear-physics corrections to this procedure — radiative corrections, weak magnetism, finite-size effects, screening — and produced the now-standard ${}^{235}$U, ${}^{239}$Pu, and ${}^{241}$Pu spectra with quoted total uncertainties around 2.5%. Mueller and collaborators in the same year extended the calculation to ${}^{238}$U, for which no ILL measurement exists, using a hybrid summation approach.

The ab-initio summation approach

The alternative is ab-initio summation: instead of using measured aggregate beta spectra, build the antineutrino prediction directly from the tabulated decay data of every individual fission fragment in nuclear-data libraries such as ENSDF or JEFF.

For each fragment, the database lists the decay Q-value, the branching ratios to different excited states, the spin and parity of parent and daughter, and the shape of the resulting beta spectrum. Multiplying by the cumulative fission yield of that fragment (which itself depends on the parent fissioning isotope) and summing over all fragments gives the predicted total antineutrino spectrum.

The advantage is conceptual cleanness: no measured beta spectrum has to be converted, every individual decay is treated explicitly, and the input data are continually updated as new nuclear-physics measurements appear. The disadvantage is systematic incompleteness: nuclear data libraries are known to be missing some short-lived high-Q-value fragments, and the listed branching ratios for many fragments are inferred from intricate analyses rather than direct measurements. The summation calculation tends to underestimate the high-energy tail of the antineutrino spectrum as a result.

Modern reactor flux models use a combination of both: conversion for the regions where the ILL measurements anchor the result, and ab-initio summation to fill in regions and isotopes where the conversion data are absent or of lower quality. Huber-Mueller is the most widely adopted such hybrid.

The reactor antineutrino anomaly

By 2011 the various reactor experiments — Bugey, Krasnoyarsk, ILL, Goesgen, Rovno, and earlier short-baseline measurements — had accumulated enough data that a global comparison with Huber-Mueller predictions could be made. The result: the world-averaged measured-to-predicted ratio sat at $0.94\pm 0.02$, a roughly 6% deficit. This is the reactor antineutrino anomaly.

Interpreted as oscillation, the deficit pointed toward a sterile-neutrino mass-squared splitting around 1-2 eV² with a sizeable mixing — coincidentally the same parameter space favoured by the gallium anomaly and the MiniBooNE excess.

But two other interpretations were also possible. The deficit could be due to a systematic underestimate of the reactor flux normalisation in Huber-Mueller, in which case there would be no sterile neutrinos — the calculation was just slightly wrong. Or the conversion approach had a hidden bias that produced a specific underestimate for some fissioning isotope.

Disentangling these required a measurement of the flux contributions from individual isotopes, not just the global average.

Daya Bay’s evolution measurement

A nuclear-reactor fuel cycle naturally provides the discriminating data. At the start of a fresh fuel cycle, the reactor burns mostly ${}^{235}$U with a small ${}^{238}$U contribution. Over the next 12-18 months, ${}^{235}$U is depleted by fission, while ${}^{239}$Pu accumulates from ${}^{238}$U neutron capture and ${}^{241}$Pu in turn builds up. By the end of the cycle the plutonium isotopes contribute roughly half the antineutrino flux.

This time-varying composition means the measured antineutrino flux and spectrum are not constant over the cycle — and the variation can be exploited to disentangle the individual isotope contributions.

Daya Bay’s evolution measurement, published in 2019, used six years of data correlated with each reactor’s instantaneous fuel composition. The analysis fit the time-varying flux and spectrum to a linear combination of per-isotope contributions, with the time-dependence of each contribution providing the lever arm.

The headline result: the measured per-fission ${}^{235}$U antineutrino yield was about 8% below the Huber-Mueller prediction, while the ${}^{239}$Pu yield was consistent with prediction (within roughly 3% statistical uncertainty). The 6% global deficit in the world average was therefore dominated by a discrepancy in the ${}^{235}$U component specifically, not by a uniform shift across all isotopes.

This is not what an oscillation explanation would predict. A sterile-neutrino-induced disappearance would affect all isotopes equally (since the disappearance is a function of energy, not of fuel composition), producing equal fractional deficits across ${}^{235}$U and ${}^{239}$Pu. The observed isotope dependence instead points to a problem in the flux modelling specifically for ${}^{235}$U — most likely a systematic error in Huber’s conversion of the ILL measurement, perhaps in the treatment of the high-energy tail or in the fissioning isotope’s particular fragment yield distribution.

The result substantially reduced the credibility of the sterile-neutrino interpretation of the reactor anomaly, though it did not by itself eliminate the possibility of a smaller sterile mixing on top of the dominant flux-modelling issue.

The 5-MeV bump

A separate and persistent feature in the measured reactor spectra is the 5 MeV bump — a 10-12% excess in the measured antineutrino flux at energies around 5 MeV relative to the Huber-Mueller prediction. The bump was first identified by RENO in 2014 and quickly confirmed by Daya Bay and Double Chooz. It is not consistent with any reasonable sterile-neutrino oscillation pattern, since it appears at all baselines.

The leading explanation is a specific mis-prediction of the beta spectra of one or more particular high-Q-value fission fragments. The conversion approach used in Huber inherently has limited sensitivity to the highest-Q-value branches because they contribute only at the high-energy end of the cumulative spectrum where statistics are poor. If the corresponding antineutrino spectra were assigned the wrong Q-values or the wrong forbidden-vs-allowed shape factors, the predicted spectrum near 5 MeV could be biased.

Daya Bay’s isotope-decomposed measurement also showed that the 5 MeV bump is dominantly in the ${}^{235}$U component, consistent with the global-deficit finding and with the interpretation that something specific about the ${}^{235}$U conversion is wrong.

What’s being done

Several efforts are improving the situation.

Direct measurements of fragment beta spectra are being made at JYFLTRAP and other isotope-separator facilities, providing fragment-by-fragment data that the summation approach can use directly. These measurements have begun to fill gaps in the nuclear-data libraries and are gradually refining the ab-initio predictions.

New ILL-style cumulative-spectrum measurements are being considered to replace the 1985 data with modern technology. These would have substantially better statistics, better resolution, and well-controlled systematic uncertainties.

Reactor experiments with on-site beta-spectrum measurement capability, such as the PROSPECT detector adjacent to the HFIR reactor at Oak Ridge, can in principle measure the antineutrino flux per fissioning isotope directly through the time-evolution and the localized geometry. PROSPECT’s ${}^{235}$U-pure HEU reactor source was particularly well-matched to this question.

Future short-baseline reactor experiments such as JUNO-TAO will provide high-statistics reference spectrum measurements with much improved energy resolution, anchoring the absolute flux prediction empirically and reducing the dependence on the model calculations.

The combined effect is that the reactor anomaly is shrinking with each refinement, and the picture is converging on the interpretation that Huber-Mueller mis-models the ${}^{235}$U component specifically while the other isotopes are fine.

Why this matters for oscillation physics

The reactor flux modelling is not just a detail. It directly affects the precision of ${\theta }_{13}$ extracted by Daya Bay, RENO, and Double Chooz — currently the most precisely measured mixing angle in the lepton sector. JUNO’s mass-ordering measurement depends critically on the absolute flux normalisation and the precise spectral shape, and a 3-5% systematic in the flux prediction translates into a substantial reduction in JUNO’s mass-ordering sensitivity. Future short-baseline searches for sterile neutrinos — at JUNO-TAO or at dedicated facilities — depend entirely on getting the absolute flux right.

The history of the reactor anomaly is also a useful cautionary tale about anomaly interpretation in general. A 6% global deficit that initially looked like new physics turned out, on closer inspection with better data, to be a 8% deficit in one specific isotope’s modelled contribution — a flux-modelling issue rather than a particle-physics signal. The sterile-neutrino space allowed by the reactor data has shrunk dramatically as the isotope-decomposition picture has clarified, leaving very little room for new physics in the reactor sector specifically.

Summary

The Huber-Mueller model is the de facto standard prediction for reactor antineutrino flux, combining Huber’s conversion of the 1985 ILL cumulative beta-spectrum measurements for ${}^{235}$U, ${}^{239}$Pu, and ${}^{241}$Pu with Mueller’s ab-initio summation for ${}^{238}$U. The quoted total uncertainty is around 2.5%, but world-average reactor data sits 6% below the prediction (the reactor antineutrino anomaly), and a $\sim$10% excess appears in the 5 MeV region. Daya Bay’s 2019 evolution measurement used the time-varying fuel composition to extract the per-isotope antineutrino yield separately, finding that the deficit is concentrated in the ${}^{235}$U component (8% below prediction) while ${}^{239}$Pu agrees with prediction. The isotope dependence is incompatible with a global sterile-neutrino oscillation and points instead to a specific flux-modelling problem for ${}^{235}$U — most likely in the conversion of the 1985 ILL measurement. Ongoing work at JYFLTRAP, PROSPECT, and the proposed JUNO-TAO will refine the flux predictions further and reduce the modelling uncertainty that limits the precision of ${\theta }_{13}$, mass-ordering, and sterile-neutrino searches in the reactor sector.