Particle physics origin of the 5 MeV bump in the reactor antineutrino spectrum?

One of the most puzzling questions in neutrino physics is the origin of the excess at 5 MeV in the reactor antineutrino spectrum. In this paper, we explore the excess via the reaction $^{13}$C$(\overline{\nu}, \overline{\nu}^\prime n)^{12}$C$^*$ in organic scintillator detectors. The de-excitation of $^{12}$C$^*$ yields a prompt $4.4$ MeV photon, while the thermalization of the product neutron causes proton recoils, which in turn yield an additional prompt energy contribution with finite width. Together, these effects can mimic an inverse beta decay event with around 5 MeV energy. We consider several non-standard neutrino interactions to produce such a process and find that the parameter space preferred by Daya Bay is disfavored by measurements of neutrino-induced deuteron disintegration and coherent elastic neutrino-nucleus scattering. While non-minimal particle physics scenarios may be viable, a nuclear physics solution to this anomaly appears more appealing.


I. INTRODUCTION
Neutrinos 1 from nuclear reactors have a prominent role in the field of neutrino physics; neutrinos were discovered using reactors as sources, and θ 13 was shown to be nonzero using detectors close to a nuclear reactor [1][2][3]. Currently, reactor experiments are grappling with two anomalous findings. The first is the ∼ 6% mismatch between the predicted and observed rates in short-baseline reactor experiments that could be explained via oscillations into additional, eV-scale neutrinos [4]. The other, more recent one is the excess of neutrino events around 5 MeV, the so-called reactor bump [5]. This spectral feature has triggered significant interest in the community in recent years [6][7][8], but has thus far eluded any quantitative explanation in terms of nuclear physics. Moreover, no attempts in terms of beyond the Standard Model (SM) physics can be found in the literature. Due to the lack of general agreement on the source of the discrepancy, it appears meaningful to explore potential for explaining the 5 MeV bump with new physics (NP).
Current reactor experiments detect neutrinos via inverse beta decay (IBD):ν e + p → n + e + . The signature used to identify these events and separate them from background is the delayed coincidence between the positron, e + -which forms the prompt signal -and the neutron capture on gadolinium and the subsequent emission of a gamma cascade of about 8 MeV total energythe delayed signal, some 50-100 µs later. The positron signal has no specific characteristics and any prompt energy deposition in the detector could be mistaken as the prompt signal. It is the delayed neutron capture signal which is essentially background free, i.e., only neutrons can cause it, which then allows to tag the prompt signal as being part of an IBD event. Under these conditions, it is clear that any new physics explanation must produce a neutron in its final state. It also has been established that the bump tracks reactor power and thus has to be caused by a particle flux from the reactor which is closely related to the number of beta decays inside the reactor. Beta decays of any known isotope extend at most to Q-values of 15-20 MeV, and thus, direct production of a neutron from any particle flux from a reactor is excluded. Therefore, a neutron which is already present in the detector material has to be liberated from its nuclear bound state.
Liquid scintillator as used in Daya Bay and NEOS is made up of carbon and hydrogen. Hydrogen has only about a 1:10,000 admixture of deuterium; thus, the only neutron source is carbon. Carbon-12 has a neutron separation energy of 18.7 MeV and thus there are no particles energetic enough coming from the reactor which can effect this reaction. On the other hand, carbon-13 has an abundance of about 1% and a neutron separation energy of only 4.9 MeV, well within the energies available in beta decay. The reaction 13 C(ν,ν n) 12 C can proceed via the ground state of 12 C or via the first excited state at 4.4 MeV, which then would de-excite immediately via the emission of a 4.4 MeV gamma, which naturally would appear as a line-like energy deposition close to the bump energy. Additional prompt energy will be supplied by the neutron kinetic energy: the neutron loses its energy by collisions with carbon and hydrogen and the recoiling nuclei will leave a scintillation signature. Due to the details of scintillation light emission, explained in more detail later, and the stochastic nature of neutron energy losses, this recoil signal will add the missing energy to match the mean bump energy as well as give the requisite width to the signal. Note, that this is a minimal scenario from a phenomenological view point and apart from the interaction causing the anomalously high 13 C scattering rate, all ingredients are well-known Standard Model physics. In particular, the bump position and width are determined by known physics.
In this paper, we will first explore the details of signal formation and demonstrate that it actually provides an excellent description of the experimental data. Then we proceed, in the spirit of low-energy effective field theory, to explore possible operators that could cause the anomalously high 13 C scattering rate and examine if other data arXiv:1803.08506v2 [hep-ph] 28 Mar 2019 rule out these operators at the required strength. We are looking for neutral-current-like operators of the form where N denotes a nucleon and L is an arbitrary Lorentz structure, i.e., L = 1, γ 5 , γ µ , γ µ γ 5 , . . .. The strength of the new interaction, , is bounded from below by the fact that the bump is about 10% of the regular IBD signal. Combined with the 1% abundance of 13 C, this implies that > 10. Invariance under SU(2) EW implies the existence of a corresponding term involving charged leptons; any such operator with the required interaction strength will be ruled out by charged-lepton data. We circumvent this restriction by introducing an additional neutrino state that is not part of an electroweak multiplet. 2 We can choose the new mixing angle and new mass splitting to be consistent with this additional neutrino being the explanation of the reactor antineutrino anomaly [4]. We argue that this constitutes a phenomenological minimal framework to seek a BSM explanation of the 5 MeV bump and we provide the quantitative demonstration that this framework indeed explains the reactor data. We also will demonstrate that any minimal incarnation of the operator in Eq. 1 is ruled out by either COHERENT [9] data or past reactor neutrino deuterium scattering results [10].

II.1. Obtaining the new physics cross section
A first-principles calculation of the 13 C(ν,ν n) 12 C ( * ) cross section for a given NP model is beyond the scope of this work. Here, we illustrate our estimation of the cross section using measurements of 13 C(e, e n) 12 C ( * ) as a proxy for the NP 13 C(ν,ν n) 12 C ( * ) cross section. In doing so, we are implicitly assuming the existence of a new vector interaction; we discuss this scenario in more detail in section IV. We assume that nuclear physics effects in the exchange of the hidden vector particle, dubbed X, are not qualitatively different from the SM photon case.
The differential cross section for 13 C(e, e n) 12 C ( * ) is parameterized in Ref. [11] by where E i is the initial-state electron energy, E n is the neutron kinetic energy, Ω f is the final-state electron solid angle, A 0 is an experimentally-determined function of the energy transferred to the nuclear system [11] and dσ Mott /dΩ f is the differential Mott cross section for an electron scattering off a heavy, point-like nucleus in QED. The Mott cross section in eq. (2) is formally divergent; we regulate it by introducing a finite mass M X for X and find We convert the (e, e n) cross section σ into the NP cross section σ for the (ν,ν n) reaction by rescaling the former by an effective coupling ζ 2 . The resulting neutron spectrum is where Eν is the energy of an incoming neutrino and dσ/dE n = 4πA 0 σ Mott . Convolving this neutron spectrum with the reactor antineutrino flux, φν, yields the cross-section-weighted neutron spectrum, φ n , where "0.006" accounts for the difference in the number of protons and carbon atoms in the linear alkylbenzene (LAB), (C 6 H 5 ) C n H 2n+1 (n ∼ 10 − 16), used by Daya Bay, as well as the natural abundance of 13 C (1.07%). Once detector effects are included, this spectrum is added to the IBD prompt-energy spectrum in order to replicate the 5 MeV bump.

II.2. Reactor Antineutrino fluxes
The 13 C(ν,ν n) 12 C * reaction requires a significant neutrino flux above Eν 9.4 MeV, though reactor experiments have measured the flux only up to 8 MeV, e.g., see Ref. [12]. The Huber-Mueller (HM) fluxes [13,14] do not extend beyond 9 MeV. On the other hand, ab initio calculations of reactor fluxes [15] predict a nonzero flux up to energies of 16 MeV. Preliminary upper bounds of the neutrino flux beyond 8 MeV have been reported by Daya Bay [16] and our fit does not exceed these bounds. We augment the HM fluxes by including an additional powerlaw component beyond 8 MeV, described by a power-law index I, such that the flux is continuous at Eν = 8 MeV.
A popular alternative explanation for the reactor bump is that the neutrino flux from 238 U has been miscalculated; there exists mounting evidence supporting a connection between the strength of the reactor bump and the 238 U fuel fraction of the reactor (see Refs. [6,17] for more details). This explanation and that proposed here are not mutually exclusive; if the flux of high-energy antineutrinos from 238 U has been underestimated, then this may provide the additional flux that we require to reproduce the bump.

II.3. Neutron thermalization
Neutrons lose energy via elastic collisions with the atoms in the scintillator; we simulate these stochastic energy losses with a simple Monte Carlo calculation. To determine the amount of prompt energy from the proton recoils, we account for quenching in the liquid scintillator, assuming nonrelativistic protons. The amount of energy produced by the scintillator is given by Birks' law [18], where k B = 0.0065 g cm −2 MeV −1 is Birks' constant for Gd-doped LAB [19] as used in Daya Bay. The quenching factor Q = E(k B )/E(k B → 0) is calculated via eq. (7) using the proton dE/dx in organic scintillator [20]. The combined effects of the stochastic proton energy distribution and quenching give rise to a finite width for the prompt energy reconstructed from this process.

III. FIT
We simulate 13 C-ν events using Monte Carlo methods for the following ranges of X masses, effective couplings and power-law indices to determine the prompt energy spectrum at Daya Bay detector(s): The smeared spectrum is compared against the ratio of measurement-to-expectation of the prompt energy spectrum reported in Ref. [6] using a standard χ 2 -test, incorporating uncertainties on and correlations between the Daya Bay data, theoretical uncertainties in the HM fluxes and experimental uncertainties on A 0 (E n ).
Considering only transitions to the excited state of 12 C, the minimum value of the χ 2 , χ 2 min , is 10.3 and occurs for {M X , log 10 ζ, I} = {1 keV, −6.075, 13}; for 21 degrees of freedom (d.o.f.), this is equivalent to p = 0.97. 3 The curve corresponding to this best-fit point is shown in fig. 1. The blue circles and blue bands correspond to the Daya Bay data and their associated uncertainties; the solid, purple line is calculated for the parameters in eq. (8) and the dotted, purple lines show the uncertainties on the neutrino spectrum and the measurements of A 0 (E n ), which have been added in quadrature. Marginalizing over the power-law index gives χ 2 min /d.o.f. = 9.1/22 (p = 0.99), compared to 35.1/24 (p = 0.07) in the absence of any NP.
When transitions to the ground state of 12 C are included, the best-fit prompt energy spectrum presents a steep upturn at low energies; the fit cannot reproduce the 5 MeV bump without overproducing low-energy events by a factor of ∼ 10, as shown in fig. 2. For the best-fit point in eq. (8), including transitions to the ground state increases χ 2 min /d.o.f. to 49.1/21 (p = 4.9 × 10 −4 ). We estimate that a O(10%) suppression of these transitions is sufficient to avoid problems. For a purely vector interaction, such a suppression is entirely ad hoc. However, if the new interaction were axial in nature, then it is likely to significantly suppress (or eliminate altogether) such transitions: the 13 C ground state has J π = 1/2 − , whereas the 12 C ground state has J π = 0 + and the first excited state in 12 C has J π = 2 + . A full calculation using the 13 C and 12 C wave functions is beyond the scope of this work.
We have also explored the effects of a threshold 4 Q 2 determine if this can lead to the required suppression of the ground state transitions, which are generally characterized by smaller momentum exchanges. The effects of a threshold can be equivalently represented by an upper bound cos θ f < 1 − Q 2 th /2E i E f . The total Mott cross section becomes ; (9) in contrast with eq. (4), in which no threshold is present. The results of this analysis are shown in Table I. The fits are obtained in a similar fashion as described for the scenario with no threshold, except that transitions to the ground state of 12 C are unsuppressed. It is clear that nonzero Q th has not appreciably lowered χ 2 min relative to the case in which Q th = 0. Together with the increased challenge of constructing a model that would exhibit this kind of behavior (since stronger NP interactions would be required), this does not seem like a viable way to suppress ground-state transitions.
Taken at face value, including the ground-state transitions still seems to be a modest improvement over the that the minimum χ 2 for a given value of M X in this range is not appreciably larger than the global minimum χ 2 . 4 The origin of such a threshold could stem for instance from an up-scattering process (which would require the presence of MeVscale partner of a sterile neutrino). min for the cases in which NP is not present; where it is present, but there is no threshold Q th in the momentum transfer; and in the cases where this threshold is 5, 10 or 15 MeV. We also show the number of degrees of freedom ("d.o.f.") in each fit, as well as the corresponding pvalue. Transitions to the ground state of 12 C are unsuppressed in these analyses. In the last row, we show the value of the χ 2 min for the case in which transitions to the ground state are totally suppressed, for comparison. SM prediction, even if ignoring them provides a better fit. However, most of the improvement to the fit in the former case comes from fitting the upturn at low energies in the Daya Bay data; the best-fit solution does little to reproduce the 5 MeV bump. We perform a similar analysis in which we omit the six lowest-energy Daya Bay data points to isolate the 5 MeV bump 5 . With this restriction, the best-fit point is identical to that of eq. (8) (independent of our previous analysis), and χ 2 min /d.o.f., having marginalized over the power-law index, is 1.
. The equivalent result in the absence of NP is 24.1/18 (p = 0.09). Importantly, this result is independent of whether or not ground-state transitions are included in the analysis. As such, the inclusion of ground-state transitions has no effect on the 5 MeV bump itself; they only cause the fit to struggle to reconcile the bump with the upturn in the low-energy data.
We briefly consider the NEOS experiment [21] which sees a similar bump at 5 MeV in their prompt energy spectrum. Including only transitions to the first excited state of 12 C and marginalizing over the power-law index, we find χ 2 min /d.o.f. = 14.2/58 (1 − p = 6.2 × 10 −10 ) at the best-fit point compared to 66.7/60 (p = 0.26) without NP. Moreover, the results of Daya Bay and NEOS are consistent with one another; analyzing both simulta- is obtained with transitions only to the first excited state of 12 C. We demonstrate the consistency between these two data sets quantitatively using a parameter goodness of fit test [22]. We calculate χ 2 PG ≡ χ 2 min,DB+NEOS − χ 2 min,DB − χ 2 min,NEOS = 0.66, which follows a χ 2 distribution with the number of d.o.f. given by 5 We emphasize that we know of no reason why these data should be considered untrustworthy. We are simply restricting our analysis to focus on the 5 MeV bump instead of the low-energy upturn.
the number of shared parameters in the fits of Daya Bay (DB) and NEOS, i.e., two. The corresponding p-value is 0.72, indicating a reasonable level of compatibility between these data sets.
In searching for a viable model for 13 C(ν, ν n) 12 C * , we first focus on the nuclear physics requirements. In order to produce a neutron and 12 C * , 13 C nuclei need to reach an excited state in scattering with neutrinos. The first two excited states of 13 C that could decay into 12 C * have either opposite parity or different spin with respect to the ground state. This precludes a scalar from being the mediator of this interaction. Nevertheless, selection rules allow for the exchange of a boson, namely a spin-2 particle in such a process. The most notable example of a spin-2 state is the graviton. However, due to the nonrenormalizability of quantum gravity, this option is not appealing. The same holds for general spin-2 particles. The simplest option is thus the exchange of a new vector particle, which we have already hinted at in section II when obtaining the 13 C(ν, ν n) 12 C ( * ) scattering cross section in eq. (4). Admittedly, we are biased towards a pure vector interaction because, as already stated, estimating the scattering cross section requires use of available data from experiments -to wit, we have used the cross section data for 13 C(e, e n) 12 C ( * ) , which is mediated by the SM photon.
We turn now to several concrete realizations of vector interactions.

IV.1. Vector interaction
We start by employing a generic dark photon coupled to the SM photon via kinetic mixing [23], where ψ is a SM fermion with charge Q ψ in units of e, F µν and F µν are the photon (A) and dark photon (X) field strength tensors, respectively. Since SM neutrinos are not charged under electromagnetism, they do not develop a coupling to the dark photon. This minimal setup cannot explain the reactor bump.
We move on to models in which SM neutrinos are charged under an extra Abelian gauge group so that they couple to X bosons. An attractive option is U (1) B−L , which is anomaly free once right-handed neutrinos are introduced [24][25][26]. This model yields ζ 2 ∼ (g B−L ) 2 2 or (g B−L ) 4 , depending on the relative sizes of the kinetic mixing and the B − L coupling constant. Unfortunately, this model faces stiff constraints, as every other SM fermion is also charged under this group; see Refs. [27,28] for details.
We are then led to an alternative Abelian symmetry, which we call U (1) X . We introduce a fourth species of neutrino (ν s ) with nonzero charge under U (1) X , and allow for nucleons to also be charged. The relevant part of the interaction Lagrangian is where g X is the U (1) X coupling constant, and Y ν , Y n are the new neutrino and neutron charge, respectively, in units of g X . When Y n = +1, U (1) X becomes gauged baryon number, U (1) B , considered in Ref. [29]. The phenomenology of an eV-scale neutrino charged under U (1) B was studied in Ref. [30]. We identify where ∆m 2 41 is the new mass-squared splitting and θ ee is the new mixing angle. For consistency with short-baseline oscillation measurements [31][32][33][34][35], we take ∆m 2 41 = 0.4 eV 2 and sin 2 2θ ee = 0.04 as our benchmark values for deriving limits and best-fits, unless otherwise stated. Given the baseline and characteristic neutrino energy, ∆m 2 41 is sufficiently large so that its oscillations average out at Daya Bay, allowing us to replace sin 2 ∆m 2 41 L 4Eν → 1 2 in eq. (12). The red (blue) band in fig. 3 represents the 68% (99%) confidence level (CL) from Daya Bay data derived by the analysis described in the previous section, after marginalizing over the power-law index I. In our analysis, we have chosen Y n = −0.65 and Y ν = 10. Constraints from pion decay (π 0 → γ + invisible) [36][37][38] imply g X 10 −2 for M X < m π and depend on neither Y n nor Y ν ; as long as the perturbativity condition Y ν g X < 4π is maintained, the new neutrino charge can be arbitrarily large. Choosing Y ν = 10 allows us to circumvent these constraints. Moreover, there exist strong bounds on this scenario from 208 Pb-n scattering [39,40]. These bounds, however, can be evaded by judiciously choosing the neutron charge to minimize the charge of 208 Pb under this interaction; the choice Y n = −0.65 allows us to do precisely this. We have also considered bounds from kaon decay [41][42][43][44] and from fifth-force searches [27,28]; these do not provide particularly stringent additional constraints.
The most relevant bounds on this model derive from neutrino scattering experiments. We show three such constraints in fig. 3: bounds from COHERENT [9] (purple), CONUS [45] (dashed violet) and deuterium scattering [10] (dot-dashed green). While COHERENT leaves a portion of the preferred Daya Bay region unconstrained, the measurement of the deuterium disintegration fully excludes this model. In the near future, this exclusion is expected to be even stronger, once the results of the CONUS experiment, which is probing coherent neutrinonucleus scattering at a reactor source, are released. In fig. 4, we show the results of a similar analysis to fig. 3, except transitions to the ground state of 12 C are included without suppression. While more of the parameter space preferred by Daya Bay lies beyond the exclusion reach of COHERENT, we do not regard this as a genuine improvement over fig. 3 for two reasons. (1) The overall quality of the fit is poorer. The minimum chi-squared, χ 2 min , is 16.8 here, whereas this value is 9.1 when transitions to the ground state are removed. This is because the fit prioritizes the low-energy upturn over the bump at 5 MeV, as mentioned above. (2) The power-law index I preferred by the fit is smaller (I ∼ 7 vs. I ∼ 13), indicating a preference for a harder antineutrino spectrum at high energies. We have not formally constrained this high-energy component in our fit, but preliminary data [16] suggest that this is inconsistent with Daya Bay observations. Therefore, we will continue to ignore transitions to the ground state of 12 C for the rest of this work.
In what follows, we briefly discuss our treatment of coherent elastic neutrino-nucleus scattering (CEνNS) and neutrino-induced deuteron disintegration.

COHERENT and CONUS
The differential cross section for CEνNS involving ν α (or ν α ) is where E ν is the neutrino energy, G F is Fermi's constant, M is the mass of target nucleus with N neutrons and Z protons, q 2 ≈ −2M E r is the momentum transfer that gives rise to nuclear recoil energy E r , F 2 (q 2 ) is the Helm form factor [46] and Q α is the species-dependent effective charge. In the Standard Model, Q α is given by where g V p = 1 2 − 2 sin 2 θ W and g V n = − 1 2 are the proton and neutron weak vector couplings, respectively, and θ W is the Weinberg angle. We ignore contributions from the axial part of the weak interactions.
The presence of a light new vector particle induces additional contributions to the effective charge that interfere with the SM. These contributions to Q α are weighted by the relevant oscillation amplitudes. Namely, where Y ν,ν is either the neutrino or antineutrino charge, as appropriate; A αα is the amplitude for ν α disappearance; and A αs is the amplitude for ν s appearance in a ν α beam. In the limit in which oscillations among the SM neutrinos are irrelevant, the oscillation amplitudes are where sin 2 θ αα ≡ |U α4 | 2 . Squaring this yields where the oscillation probabilities are given by We emphasize that the only nonzero new mixing parameter that we consider is sin 2 2θ ee , i.e., muon neutrinos do not mix with the fourth flavor.
We wish to determine the exclusion reach of CO-HERENT and the sensitivity of CONUS. We simulate COHERENT following the procedure described in Refs. [47,48]; for CONUS, we follow Ref. [48]. For most values of Y n , COHERENT excludes most of Daya Bay's preferred parameter space. However, choosing the neutron charge to minimize the U (1) X charges of the constituent nuclei allows for the bound to be loosened by a factor of ∼ 10. Since the relevant isotopes for COHER-ENT are 133 Cs and 127 I, Y n = −0.65 -which we had previously chosen to evade 208 Pb-n scattering constraintsprovides modest relief (see eq. (12)).
We show the 99% CL exclusion limit for COHER-ENT in purple in fig. 3. If we had selected the best-fit point from Ref. [35], (∆m 2 41 , sin 2 2θ ee ) ∼ (1.3 eV 2 , 0.04), then COHERENT would have ruled out the entire region preferred by Daya Bay. Given COHERENT's characteristic neutrino energy (∼ O(30) MeV) and baseline (19.3 m), the event rate approximately scales as g 4 X sin 2 2θ ee (∆m 2 41 ) 2 for ∆m 2 41 0.5 eV 2 . The COHER-ENT constraint is weakened as ∆m 2 41 is decreased, but our willingness to do so in order to avoid this constraint is dictated by reactor data [35]. Our benchmark parameters -(∆m 2 41 , sin 2 2θ ee ) ∼ (0.4 eV 2 , 0.04) -are allowed at 95% CL by Ref. [35] and permit us to marginally evade this constraint.
While the sensitivity of COHERENT is slightly diminished for Y n = −0.65, the CONUS detector is made of natural germanium, for which the proton-to-neutron ratio is ∼ 0.8. Consequently, each of the five naturallyoccurring germanium isotopes has nonzero charge under U (1) X , so it is not possible to simultaneously provide relief from both Pb-n scattering and COHERENT constraints while also limiting the sensitivity of CONUS. The 99% CL sensitivity of CONUS is shown in dashed violet in fig. 3.

Neutrino-induced deuterium disintegration
With the COHERENT limit not fully excluding the preferred Daya Bay region, and with CONUS not yet having released any results, only the bound from the neutrino-induced deuterium disintegration currently rules out the exchange of a vector boson. We briefly discuss how this limit has been calculated.
We have calculated the contributions of our NP model to the cross section for inelastic neutrino-deuterium scattering at reactors, following the calculations of Refs. [49,50]. In appendix A, we outline the computation of the cross section in the SM and demonstrate how the effects from NP can be included.
To obtain the limit, we first calculate the flux-weighted cross section at 11.2 m from a nuclear reactor, as in Ref. [10], for our benchmark fourth neutrino parameters using the HM fluxes for a range of g X and M X . The results are then compared against the same for the purely SM case. Ref. [10] reports a measurement of the ν scattering rate consistent with the SM value at the ∼ 10% level; this allows us to exclude parts of the {g X , M X } parameter space. In fig. 3, we show in light green the contour along whichσ SM+NP = 2σ SM in the {g X , M X } plane, whereσ SM (σ SM+NP ) is the SM (SM+NP) fluxweighted cross section; this is intended to be a conservative limit, given the complexity of the underlying measurement. This result clearly excludes the vector-bosonexchange explanation of the reactor antineutrino excess.

IV.2. Axial interaction
Next, we consider a spin-1 boson X with purely axial couplings to the fourth neutrinos and to nucleons, L ⊃ g X X µ (Y ν ν s γ µ γ 5 ν s + pγ µ γ 5 p + Y n nγ µ γ 5 n) . (20) Integrating out X leads to the effective Lagrangian of the form where we have selected the left-handed part of the neutrino current because any neutrino present in an experiment will arise from a neutrino (antineutrino) produced in weak-interaction processes and will thus be initially left-handed (right-handed), to leading order.
Naively, an axial-vector interaction may be able to relax constraints from COHERENT and CONUS; for similar vector and axial couplings to nucleons, the axial contribution to the cross section is roughly a factor of 1/A smaller than its vector counterpart [51,52]. Hence, if the NP vector couplings were to vanish, then the axial interaction may be significantly less constrained than a comparably-sized vector interaction. We have performed such an analysis and have found a factor of few relaxation of the g X limit. This qualitatively does not alter the picture given in section IV.1. It is interesting to note that it is possible to make either the COHERENT (CONUS) bound vanish, at leading order, by demanding the proton (neutron) to be uncharged under this new axial interaction.
Despite the more positive outlook at CEνNS experiments compared to the pure vector interaction, the bound from deuterium disintegration turns out to be stronger than its vector-interaction counterpart. We do not review here the full derivation, but point the interested reader to appendix A, where the treatment for the vector interaction is outlined. We emphasize that only the coefficients representing the axial-vector isoscalar and isovector couplings (C (0) A ) get modified in this case. As mentioned previously, a full calculation of the 13 C(ν, ν n) 12 C ( * ) cross section for an axial-vector mediator is well beyond the scope of this work. However, it is likely that neutrino-induced deuterium disintegration would place a cripplingly strong bound on such a scenario.

IV.3. Magnetic interaction
Lastly, we consider a magnetic-magnetic interaction between the new neutrinos and nuclei. We introduce an effective interaction of the form where µ ν and µ N are the neutrino and nucleon magnetic dipole moments, respectively. We parameterize the latter as µ N = g X µ SM N , where µ SM N is the SM magnetic dipole moment with the fundamental charge e factored out (µ SM p = 2.79/2m N and µ SM n = −1.91/2m N ). The neutrino magnetic moment is defined to be µ ν ≡ g X Y ν /M X .
To estimate the constraint from Daya Bay, we follow the example of Ref. [53], wherein the νd → νnp cross section is calculated for nonzero neutrino magnetic moments using the γd → np cross section to sidestep complications arising from nuclear structure. Of course, the 13 C(X, n) 12 C ( * ) cross section is not measured; we use the 13 C(γ, n) 12 C ( * ) cross section instead [54] modified by a factor of 2 3 (g X /e) 2 , to account for (1) the difference in coupling strengths, and (2) the difference in the number of allowed spin states of our massive vector boson. This forces us into assuming that the longitudinal mode of our massive boson would not contribute appreciably to 13 C(X, n) 12 C ( * ) . The experimentally determined total 13 C(γ, n) 12 C ( * ) cross section does not reveal how much of this cross section comes from M1 scattering and how much from E1 scattering. We assume that M1 scattering dominates; if this is untrue, then the Daya Bay data will prefer larger values of g X than this optimistic analysis.
With these assumptions and by following Ref. [53], we obtain the following differential cross section: where E (E ) is the initial (final) antineutrino energy and σ γ is the 13 C(γ, n) 12 C ( * ) cross section. Note that the final neutrino energy satisfies where S n is the neutron removal energy of 13 C, E ex is the excitation energy of the final 12 C (either 0 or ∼ 4.4 MeV), and E nuc is the kinetic energy carried off by the 12 C-neutron system. We assume that the neutron carries off most of this energy and, therefore, that eq. (23) is also the neutron recoil spectrum.
The red (blue) band in fig. 5 is the 68% (99%) confidence level preferred by the Daya Bay data. We include only transitions to the first excited state of 12  For the CEνNS cross section, we extend the formalism for nonrelativistic dark matter-nucleon effective field theory presented in Refs. [55,56] to account for the fact that neutrinos are ultrarelativistic; a detailed calculation will be presented in Ref. [57]. We have verified that we recover the usual CEνNS cross section from this framework in the absence of NP. The result is similar to the one pre-sented in section IV.1 for a vector interaction. In fig. 5, we see that COHERENT (purple) does little to challenge the parameter space preferred by Daya Bay. We have shown the result for our benchmark fourth neutrino parameters, in the interest of consistency, but have determined that even for the best-fit point of Ref. [35], having a larger ∆m 2 41 ( = 1.3 eV 2 ), COHERENT provides only a modest exclusion of the high-mass (M X 10 MeV) portion of the Daya Bay preferred region. However, CONUS (dashed violet) will have the power to fully exclude this model.
For deuterium disintegration, we repurpose the results of Ref. [58] to determine the cross section. We use a modified version of the M1 matrix element given in their eq. (11), to account for differences in normalization convention and the effects of a nonzero M X : where n 1 (n 2 ) is the unit vector along the direction of the incident (outgoing) neutrino, which respectively have energies E and E f . Ref. [58] only includes the isovector contributions to the scattering cross section because these would ultimately dominate. Additionally, |M nuc | 2 incorporates the effects of the structure of deuterium; it is given by Eq. (12) of Ref. [58] and depends on E k , the kinetic energy of the relative motion of the final-state proton and neutron.
The differential cross section reads where µ SM ≡ 2m N (µ SM p − µ SM n ) is the isovector nucleon magnetic moment, m N is the nucleon mass and B is deuteron binding energy. The function f (z) reads where we define z ≡ M 2 X /(2EE f ). From eq. (25), it is straightforward to calculate the flux-weighted cross section and obtain the limit in a similar fashion as described at the end of section IV.1. From fig. 5, we observe that the deuterium disintegration limit (dot-dashed green) again rules out the entirety of the parameter space preferred by Daya Bay. Hence, this model is strongly disfavored.

V. SUMMARY AND CONCLUSIONS
In this work, we have studied whether the 5 MeV neutrino excess at reactors can be addressed with physics beyond the SM -more precisely, via the interplay between nuclear physics and new particle physics. By injecting at least 9.4 MeV of energy into the nucleus, 13 C may transition to the lowest-lying excited state of 12 C which subsequently de-excites. The resulting 4.4 MeV photon, together with energy deposited during neutron thermalization, can successfully explain the excess. We introduced additional, eV-scale neutrinos and assumed them to have hidden interactions with nucleons to allow them to scatter off of 13 C nuclei.
We considered three minimal types of the interactions: vector, axial and magnetic. In all cases, the parameter space preferred by Daya Bay is excluded by neutrinoinduced deuterium disintegration, is in severe tension with the data collected by COHERENT and will be strongly challenged by CONUS data in the near future. While non-minimal NP scenarios could explain the reactor bump, this work suggests that the antineutrino excess in the reactor data is more likely to be of nuclear physics origin.
Interestingly, is possible to probe whether or not 13 C(ν, ν n) 12 C * is responsible for the reactor bump in a model-independent way. This reaction is available at all organic scintillator detectors (liquid or solid) and the prompt signal is a combination of proton recoils (about 20%) and the 4.4 MeV de-excitation photon. In singlevolume detectors like Daya Bay and NEOS, the 4.4 MeV photon is well contained; neither employs pulse shape discrimination, so these events pass as genuine IBD events.
PROSPECT [59][60][61] and STEREO [62][63][64] both rely on pulse shape discrimination to manage backgrounds, so the small proton recoil contribution will likely cause these events to be rejected as background. Also, the 4.4 MeV photon would deposit energy across many detector cells; it seems likely that no bump would be observed in the standard analysis. The relatively long distance over which the 4.4 MeV photon deposits its energy, compared to a positron in an IBD event, is likely to have these events rejected in finely segmented detectors like SoLid [65] and DANSS [66,67]. Curiously, DANSS seems not to observe a bump 6 , in accordance with this prediction. A dedicated analysis in these segmented detectors could likely identify this process 7 .
The analogous reaction 17 O(ν,ν n) 16 O * is unlikely to be observed at water Cerenkov detectors. First, the natural abundance of 17 O (0.03%) is significantly lower than that of 13 C (1.07%). Second, while the neutron separation energy of 17 O (4143.1 keV) is slightly less than that of 13 C, the lowest-lying excited state of 16 O has an excitation energy of 6049.6 keV [70]; populating this state would require a significant antineutrino flux beyond E ν 10 MeV, and we have seen that this flux dies off quite quickly. Even though this could, in principle, give rise to an analogous bump in the prompt energy spectrum at ∼ 6 MeV in water Cerenkov detectors, it would be substantially more difficult to detect. Therefore, neither Gd-doped Super-K [71] nor WATCHMAN [72] should see a bump.
6 DANSS collaboration recently announced the observation of an excess around 5 MeV with a low statistical significance. They also pointed out that more effort in the direction of calibration is required in order to robustly claim the excess (see [68] for details). 7 Since this work first appeared on arXiv, Ref. [69] reported an observation of a bump at 5 MeV in the segmented Gösgen experiment.
where E f and Ω f are the final-state neutrino energy and solid angle, respectively. The functions W LO 0 and W LO 1 are hadronic structure functions, calculated in Refs [49,50]. For completeness, we list the functions that enter into this cross section: 4M N γp π(p 2 + γ 2 ) 2 1 − |q| 2 (p 2 + 3γ 2 ) 6(p 2 + γ 2 ) 2 , (A5) In these formulae, M N ≈ 938 MeV is the nucleon mass, B = 2.2245 MeV is the deuteron binding energy, a ( 1 S0) = −23.7 fm is the deuteron scattering length, C where θ W is the Weinberg angle, ∆s is the strange-quark contribution to the proton spin and g A is the πN N coupling constant.
The new vector interaction we introduce induces additional contributions to C (I=0,1) V . Starting from the effective Lagrangian we determine the NP isoscalar and isovector couplings to be , (A20) where we have used that q 2 = −4E i E f sin 2 θ f 2 . Note that a negative sign has been absorbed via Y ν = −Y ν . We incorporate the SM and NP by adding C where we need only to consider oscillations involving electron (anti)neutrinos. Squaring this and using the expressions in eqs. (16) and (17) To calculate the total cross section, eq. (A1) must be integrated over the region