Nuclear de-excitations in low-energy charged-current $\nu_e$ scattering on $^{40}$Ar

Background: Large argon-based neutrino detectors, such as those planned for the Deep Underground Neutrino Experiment (DUNE), have the potential to provide unique sensitivity to low-energy ($\sim$10 MeV) electron neutrinos produced by core-collapse supernovae. Despite their importance for neutrino energy reconstruction, nuclear de-excitations following charged-current $\nu_e$ absorption on $^{40}$Ar have never been studied in detail at supernova energies. Purpose: I develop a model of nuclear de-excitations that occur following the $^{40}\mathrm{Ar}(\nu_e,e^{-})^{40}\mathrm{K}^*$ reaction. This model is applied to the calculation of exclusive cross sections. Methods: A simple expression for the inclusive differential cross section is derived under the allowed approximation. Nuclear de-excitations are described using a combination of measured $\gamma$-ray decay schemes and the Hauser-Feshbach statistical model. All calculations are carried out using a novel Monte Carlo event generator called MARLEY (Model of Argon Reaction Low Energy Yields). Results: Various total and differential cross sections are presented. Two de-excitation modes, one involving only $\gamma$-rays and the other including single neutron emission, are found to be dominant at few tens-of-MeV energies. Conclusions: Nuclear de-excitations have a strong impact on the achievable energy resolution for supernova $\nu_e$ detection in liquid argon. Tagging events involving neutron emission, though difficult, could substantially improve energy reconstruction. Given a suitable calculation of the inclusive cross section, the MARLEY nuclear de-excitation model may readily be applied to other scattering processes.


I. INTRODUCTION
Core-collapse supernovae are exceptionally intense sources of tens-of-MeV neutrinos and antineutrinos of all flavors. In a typical supernova, about 10 58 neutrinos are released in a burst lasting tens of seconds. Although the first observation of supernova neutrinos by the Kamiokande-II [1], Baksan [2], and IMB [3] detectors in 1987 yielded a total of only two dozen events, the scientific impact of this dataset has been tremendous, leading to numerous publications on a wide variety of subjects [4][5][6][7][8][9]. In the years since first detection, a worldwide network of large neutrino experiments, most built primarily for other applications, stands ready to perform a second, high-statistics measurement if a core-collapse supernova should occur within our galaxy [10]. Due to the slow rate of galactic core-collapse supernovae (estimated to be about 1.6 per century [11]), the prospect of such a measurement represents a rare but valuable opportunity to shed light on the details of core-collapse and nucleosynthesis, study neutrinos under extreme conditions, search for evidence of physics beyond the Standard Model, and explore many other topics [12][13][14][15].
A full realization of the physics potential of the next galactic core-collapse supernova will require a simultaneous measurement of neutrinos of all flavors. While detectors based on water and hydrocarbon scintillator will primarily detectν e via inverse beta decay (IBD) ν e + p → e + + n , (1) * gardiner@fnal.gov liquid-argon-based detectors are anticipated to provide unique sensitivity [16,17] to ν e via the charged-current (CC) reaction ν e + 40 Ar → e − + 40 K * which dominates the expected signal at supernova energies.
Within the decade, the Deep Underground Neutrino Experiment (DUNE) will begin operating four ten-kiloton liquid argon time projection chambers (LArTPCs) with the primary goals of studying longbaseline oscillations of accelerator neutrinos, searching for proton decay, and measuring the ν e flux from a galactic core-collapse supernova if one should occur during the lifetime of the experiment [18]. Initial studies of the sensitivity of DUNE to supernova neutrinos, performed by the collaboration itself [16] and by smaller groups (e.g., ref. [19]) show promise, and the potential exists for measurements by DUNE of other low-energy neutrinos, notably those produced by the Sun [20].
In addition to DUNE, three sub-kiloton LArTPCs, SBND [21,22], MicroBooNE [23], and ICARUS [24], are currently operating or being installed in the Booster Neutrino Beam at Fermilab. A joint effort between the three experimental collaborations, known as the Short Baseline Neutrino (SBN) program [25,26], will perform precision measurements of neutrino oscillations. In addition to this primary mission, the SBN detectors will pursue a variety of other physics measurements and are expected to be sensitive to supernova neutrinos. To ensure that data from a core-collapse supernova would be fully recorded over the ∼10 s burst, the MicroBooNE collaboration op-erates a first-of-its-kind continuous readout stream and has demonstrated its capabilities via reconstruction of Michel electrons produced by decays of cosmic muons [27].
While much remains to be done to fully exploit the lowenergy capabilities of LArTPCs, a first demonstration by the ArgoNeuT [28] experiment of reconstruction of MeV-scale activity due to accelerator-neutrino-induced neutrons and de-excitation γ-rays achieved a detection threshold of ∼200-300 keV [29]. These encouraging initial results have prompted further experimental work by MicroBooNE [30] and multiple simulation-based studies considering the implications for reconstruction of both high-and low-energy physics events [31,32].
In future analyses of supernova neutrino data, the event-by-event reconstructed neutrino energy obtained by each detector will be of primary interest. For IBD events in water or scintillator, because only a single hadronic final state (a free neutron) is accessible at tensof-MeV energies, a measurement of the outgoing positron energy is sufficient to reconstruct the antineutrino energy with high accuracy. Due to the use of a complex nuclear target (argon) in LArTPCs, however, various nuclear transitions may occur in response to CC ν e absorption, and thus a one-to-one mapping (up to nuclear recoil) between the neutrino and electron energies no longer exists.
To fully reconstruct the neutrino energy in the argon case, the reaction Q-value, i.e., the energy imparted to the nuclear transition, must be inferred by detecting the nuclear de-excitation products. For transitions to bound nuclear energy levels, the neutrino energy is in principle fully recoverable by measuring the energies of all deexcitation γ-rays in addition to the primary electron. For transitions to unbound nuclear states, a model is needed to correct for missing energy associated with undetected nuclear fragments. In practice, an experimental analysis that attempts to isolate the simpler bound transitions will also need a detailed de-excitation model in order to estimate the purity of the event selection.
Experimental data have not yet been obtained for neutrino-argon cross sections in the supernova energy regime (and only a few measurements with limited precision are available for any nuclear target [33, table 3]). Despite this, a substantial literature exists for theoretical calculations of the 40 Ar(ν e , e − ) 40 K * process. A review through 2018 is provided in ref. [34, sec. 7.1], with a notable recent addition being two publications [35,36] which employ a Continuum Random Phase Approximation (CRPA) model to study this cross section above the nucleon emission threshold. While highly useful for providing competing estimates of event rates in DUNE and other argon-based detectors, all published calculations for this cross section to date share the limitation of being fully inclusive, i.e., predictions are made that consider only the kinematics of the outgoing electron. At very low neutrino energies, where only transitions to bound nuclear states are pos-sible, this is not problematic: measured de-excitation γray branching ratios exist for many levels of the daughter 40 K nucleus, and missing data may be addressed using straightforward theoretical estimates. However, above ∼15-20 MeV, kinematic access to unbound nuclear states becomes appreciable, and a detailed treatment of the competition between various de-excitation channels (including emission of both γ-rays and nuclear fragments) is needed.
Although such a treatment has not previously been provided for 40 Ar, detailed modeling of nuclear deexcitations induced by low-energy neutrino interactions has been pursued for a number of other nuclei [37][38][39][40][41][42]. A universal assumption made by all of these calculations (as well as the present work) is that of compound nucleus formation: following the primary interaction, the nucleus is left in a thermally-equilibrated excited state that decays independently of the details of its formation process. While further work is needed to fully investigate the adequacy of this assumption for low-energy neutrinonucleus reactions, both theoretical calculations [43,44] and electron scattering data [45] suggest that compound processes dominate over the direct nucleon knock-out important at higher energies.
In this paper, I present the first calculations at tens-of-MeV energies for cross sections for exclusive final states of the reaction 40 Ar(ν e , e − ) 40 K * , emphasizing the role of nuclear de-excitation processes. In section II, I develop a simple model for the inclusive differential cross section, relying on approximations that work best at low momentum transfers. The derivation in section II fully determines the cross section up to the values of two nuclear matrix elements, B(F) and B(GT), which are considered in section III. While relevant neutrino scattering data are currently unavailable, the needed values of these matrix elements at low excitation energies may be extracted from measurements of related processes. I supplement these measurements with the results of a theoretical calculation at high excitation energies to obtain a full treatment of the inclusive cross section. In section IV, I describe a detailed model of nuclear de-excitations that can be used together with the inclusive cross section to obtain predictions for exclusive final states. In section V, I present example results calculated using the models developed in the previous sections.
To enable practical calculations that have already helped to inform studies of DUNE's sensitivity to supernova neutrinos [16], all of the physics models described herein have been implemented in a new Monte Carlo event generator called MARLEY (Model of Argon Reaction Low Energy Yields). All results shown in this work may be reproduced using version 1.2.0 of MARLEY [46], which is publicly available as an open-source software project [47]. Documentation of the technical details of MARLEY and usage instructions are available in ref. [48].
Due to the compound nucleus assumption, the MARLEY de-excitation model may easily be used in the future to provide exclusive predictions for a more refined calcula-tion of the inclusive CC ν e absorption cross section for 40 Ar. A similar approach to modeling de-excitations for other reaction modes and target nuclei is likewise possible, and I welcome potential collaboration on this topic. Prospects for improving MARLEY predictions beyond the proof-of-concept reported here are briefly considered in section VI.

II. INCLUSIVE CROSS SECTION MODEL
For momentum transfers that are small compared to the W boson mass, the tree-level amplitude M for inclusive charged-current neutrino-nucleus scattering may be represented diagramatically as . (3) The corresponding differential cross section may be written in the form where q = k − k = p − p is the four-momentum transfer, Q 2 ≡ −q 2 , G F is the Fermi constant, V ud is the Cabibbo-Kobayashi-Maskawa matrix element connecting the up and down quarks, Mandelstam s is the square of the total energy in the center-of-momentum (CM) frame, and m i is the mass of the initial-state nucleus. Discussion of the Coulomb correction factor F C is deferred to section II B. The leptonic (L µν ) and hadronic (W µν ) tensors are defined by and Here m is the mass of the final-state lepton, J i (J f ) is the initial (final) nuclear spin, and M i (M f ) is the third component of the nuclear spin in the initial (final) state. Under the impulse approximation, the nuclear matrix element may be written in coordinate space as where q is the three-momentum transfer and the sum runs over all A nucleons. The weak current operator j µ (n) is understood to act only on the nth nucleon, as is the position operator x(n). The state vectors in eq. (8) are normalized relativistically, i.e., where E i (E f ) is the total energy of the nucleus in the initial (final) state. Equation (4) contains an implied sum over the accessible final nuclear states.

A. Allowed approximation
The full expression for the single-nucleon weak current operator j µ is well-known and is given in ref. [34] among other places. For this study, however, I evaluate the current operator in the allowed approximation, which combines the long-wavelength (q → 0) limit and the slow-nucleon limit (in which the momentum of the initial struck nucleon is neglected compared to its mass).
Under this approximation, the weak current operator reduces to the simple form where j 0 is the time component and the three Cartesian spatial components are labeled with a ∈ {1, 2, 3}. The time component of the nuclear matrix element N µ is given by while the spatial components may be written in spherical coordinates as is the vector (axialvector) weak coupling constant of the nucleon. The Fermi and Gamow-Teller (GT) operators are defined by where σ is the Pauli vector, and t − , the isospin-lowering operator, converts a neutron into a proton. Double bars ( ) denote matrix elements which have been reduced via the Wigner-Eckhart theorem. Equations (11) and (12) may be used to evaluate the elements of the hadronic tensor W µν . Under the allowed approximation, these become where the reduced Fermi and Gamow-Teller matrix elements are given by The state vectors labeled using the nuclear spin (J i or J f ) have unit norm. The reduced matrix elements satisfy the spin-parity selection rules and Combining the results above leads to the following expression for the allowed approximation differential cross section in the CM frame: Here E , p , θ , and β = E /|p | are, respectively, the total energy, 3-momentum, scattering angle, and speed of the final-state lepton. The factor E i E f /s accounts for nuclear recoil and is commonly neglected.
In the CM frame, the particle energies are independent of the scattering angle θ . As a result, integration of the total cross section is trivial and leads to the expression (23) As was the case for eq. (4), the cross section formulae in eqs. (22) and (23) contain an implicit sum over nuclear final states.

B. Coulomb corrections
Final-state interactions (FSIs) of the outgoing charged lepton with the Coulomb field of the nucleus have a significant effect on the cross section at low energies. While a detailed treatment of Coulomb FSIs is achievable via the distorted-wave Born approximation, a much more convenient approximation scheme based on the work of Engel [49] is typically used, e.g., in refs. [35,50,51].
Under this approach, the Coulomb correction factor F C that appears in eqs. (4), (22) and (23) is calculated using either the Fermi function [52,53] or the modified effective momentum approximation (MEMA) [49]. Since the former is known to overestimate Coulomb corrections at high lepton energies while the latter does the same at low energies, the smaller of the two alternatives is always chosen. This amounts to defining the Coulomb correction factor as where the Fermi function is given by and In eq. (25) the quantity S is defined in terms of the fine structure constant α by where Z f is the proton number of the final nucleus. The nuclear radius (in natural units) may be estimated as and the Sommerfeld parameter η is given by where z is the electric charge (in units of the elementary charge) of the final-state lepton. Typical presentations of the correction factors defined in eqs. (25) and (26) neglect the small recoil kinetic energy of the final nucleus in the laboratory frame. This allows the use of expressions for F Fermi and F MEMA which are derived in the rest frame of the final nucleus. I opt instead for Lorentz-invariant forms of the correction factors written in terms of the relative speed β rel of the two final-state particles: [54] The symbols E and K from eq. (26) denote, respectively, the total energy and momentum of the outgoing lepton in the rest frame of recoiling nucleus: The effective values of these variables are those that exist in the presence of the nuclear Coulomb potential which is approximated by that at the center of a uniformly-charged sphere: It should be noted that, as originally presented [49], the MEMA also involves modifying the value of the momentum transfer used to evaluate the amplitude M. Since the cross section treatment presented here involves use of the long-wavelength limit q → 0, however, I neglect this additional correction.

III. ALLOWED NUCLEAR MATRIX ELEMENTS
Despite sustained community interest and a concrete proposal by the CAPTAIN experiment [55] to perform a direct measurement, no experimental data for tens-of-MeV neutrino scattering on argon are currently available. However, in recent decades, three separate experiments have obtained values of the allowed matrix elements B(F) and B(GT) by considering related physics processes.
The first two experiments, performed in the late 1990s by separate teams working at GSI [56] and GANIL [57], both sought to study CC ν e absorption on 40 Ar by mea-suring beta decays of its mirror nucleus 40 Ti: In the limit of perfect isospin symmetry, the matrix element describing a 40 Ti beta decay transition to a specific 40 Sc level is equal to the matrix element accessing the level's isobaric analog in 40 K via CC ν e scattering on 40 Ar. The main difficulties in applying this technique to neutrino cross sections are (1) the beta decay Q-value limits the maximum excitation energy that may be studied, and (2) energy levels in the beta decay daughter nucleus ( 40 Sc) must be matched to their analogs in the final-state nucleus for neutrino scattering ( 40 K). The third experiment [58], performed about a decade later at IUCF, extracted B(GT) values from measurements of (p,n) scattering on 40 Ar. The extraction technique relied on the observation, first put forward in 1980 [59] and subsequently refined [60][61][62], that the (p,n) cross section at very forward angles (θ ≈ 0 • ) for ∼100 MeV protons is approximately proportional to the allowed matrix elements B(F) and B(GT). While subject to some unique difficulties of its own (see, e.g., ref. [63, sec. 4.2]), this approach overcomes key limitations of 40 Ti beta decay: transitions to excited levels of 40 K may be studied directly at energies higher than the mirror beta decay Q-value.

A. Re-evaluation of existing measurements
Reasonable attempts were made in the original publications describing these measurements to assign the extracted matrix elements to known 40 K levels satisfying the spin-parity selection rules in eqs. (20)- (21). That is, the 40 K isobaric analog state accessed via a Fermi transition must have J π = 0 + , while GT transitions may only populate levels with J π = 1 + . However, in light of new 40 K level data that became available in 2017 [64], I revisited the level assignments for all three measurements.
The results of this re-evaluation are shown in table I. Level energies (keV) and spin-parity assignments retrieved from the ENSDF database [65] are listed in the first and second columns, respectively. Excitation energies (for either 40 Sc or 40 K as appropriate) and matrix element values are listed in the following columns for each of the three experimental measurements. In the case of the (p,n) scattering data, the matrix element values provided in the original paper [58] have been scaled by a factor of g 2 A = 1.26 2 . This scaling was done because the definition of B(GT) used by the experiment does not include the axial-vector weak coupling constant g A . The specific value g A = 1.26 was chosen for consistency with the one assumed in the experimental analysis. Figure 1 shows the Gamow-Teller strengths obtained by the two 40 Ti beta decay experiments mentioned previously. Excitation energies of analog states in 40 K, represented on the horizontal axis, are chosen to match the 2289.868 (11) 1 + 2287 (10)  a Parenthesized values are based upon weak arguments [66]. b The data tabulated in ref. [58] were multiplied by g 2 A = 1.26 2 to obtain the B(GT) values shown here. c A nuclear level with parity Π and spin J has natural parity if Π = (−1) J . Otherwise it has unnatural parity. d Another candidate 40 K level for this transition has Ex = 3663.88 keV and J π = (1 − , 2, 3, 4 + ). e This level is the isobaric analog of the 40 Ar ground state. f Gamow-Teller transitions are assumed for all levels other than the isobaric analog state.  40 Ti β + Bhattacharya et al. [57] (p,n) Bhattacharya et al. [58] FIG. 2: Comparison of the Gamow-Teller strengths B(GT) measured using 40 Ti beta decay (see fig. 1) with those obtained using a measurement of 0 • (p,n) scattering.
assignments made in table I. The vertical axis is inverted for the second dataset to facilitate comparisons. Rough consistency is seen between the two measurements, although the results reported by Liu et al. involve several more nuclear levels. Figure 2 uses a similar format to individually compare each beta decay measurement to the GT strengths extracted from (p,n) scattering. Substantially more fragmentation of the strength is seen in the beta decay data, and there are areas of significant tension. For instance, the two experimental methods disagree on whether the GT strength to the 40 K level at 2.3 MeV is larger or smaller than the strength to the level at 2.7 MeV.
Differences between the beta decay and (p,n) data were examined in detail by Karakoç et al. in 2014 [67]. Based on a combination of theoretical calculations and a currently unpublished 40 Ar(h, t) 40 K measurement, the authors argued that the (p,n) data should be preferred over the 40 Ti beta decay data for calculations of CC ν e absorption cross sections on argon. Rather than attempt to adjudicate between the conflicting datasets, I have opted to allow each of the three measurements to be used as a source of B(GT) values in MARLEY cross section calculations.

B. Extension to higher excitation energies
Beyond the maximum excitation energy of ∼8 MeV probed by the experiments mentioned in the previous section, the presence of additional Gamow-Teller strength is predicted by the model-independent Ikeda sum rule [68]. This rule states that the summed GT strength S − GT (S + GT ) over all possible nuclear final states for CC ν e (ν e ) absorption satisfies the relation where N i = 22 (Z i = 18) is the neutron (proton) number of the initial 40 Ar nucleus. Equation (35) implies a minimum possible value of S − GT = 12 g 2 A ≈ 19 for the integrated GT strength associated with the reaction 40 Ar(ν e , e − ) 40 K * . Comparing this value to the measured total GT strengths listed in the final row of table I reveals that the majority of the expected GT strength for CC ν e absorption on 40 Ar is unmeasured and associated with transitions to high-lying, nucleon-unbound states of 40 K.
To supplement the experimental measurements with an estimate of the remainder of the GT strength, I rely on a Quasiparticle Random Phase Approximation (QRPA) calculation by Cheoun, Ha, and Kajino [69]. To avoid double-counting, theoretical GT matrix elements are included with the experimental ones only for excitation energies at which the integrated QRPA GT strength exceeds the experimental total.  In addition to extracting GT strengths, both beta decay experiments measured values of the Fermi matrix element B(F). Under the approximation that isospin is a good quantum number, this matrix element is expected to have the value and to connect the ground state of 40 Ar to a single 0 + isobaric analog state in 40 K, which has been identified as the level with excitation energy E x = 4.3837 MeV. Since the experimental data are fully consistent with these expectations, I adopt the value B(F) = 4 for this transition in all three sets of MARLEY matrix elements. The known Fermi transition is represented in each panel of Figure 3 by a white bar with horizontal hatch marks. Transitions to all other nuclear levels are assumed to proceed via the Gamow-Teller operator. Although the differences become important at neutrino energies near threshold, the observables predicted in this paper are largely insensitive to the choice between the three sets of MARLEY matrix elements for the 40 Ar(ν e , e − ) 40 K * reaction. For definiteness, all MARLEY calculations shown in this work (see section V) are obtained using the ve40ArCC Bhattacharya1998.react input file, which contains the matrix elements shown in the middle panel of fig. 3.

IV. NUCLEAR DE-EXCITATION MODEL
To model nuclear de-excitations following CC ν e absorption on 40 Ar, I rely on the observation that, due to the selection rules in eqs. (20) and (21), each nuclear final state accessed by the neutrino interaction has a well-defined excitation energy, spin, and parity. Distinct treatments are used for bound and unbound nuclear states, with the latter being those for which the excitation energy exceeds the separation energy for at least one nuclear fragment with mass number A ≤ 4. Separation energies are computed in MARLEY using atomic and particle mass data from refs. [70,71]. Untabulated atomic masses are estimated using the liquid drop model of Myers and Swiatecki [72,73]. The full de-excitation cascade is treated as a sequence of binary decays while neglecting the possibility of fission and emission of heavy nucleon clusters (A ≥ 5).

A. Bound states: discrete level data
De-excitations of bound nuclear states are handled in MARLEY using a set of nuclear structure data files originally prepared for use with version 1.6 of the TALYS nuclear reaction code [74][75][76] [77]. These in turn are based on the level schemes included in version 3 of the Reference Input Parameter Library (RIPL-3) [78]. For a large number of nuclei, the files provide tables of discrete nuclear energy levels including their excitation energies, spin-parities, and de-excitation γ-ray branching ratios. Missing experimental measurements of these quantities are supplemented by theoretical estimates. Although in-ternal conversion coefficients are provided in the TALYS data files, this process is neglected in MARLEY 1.2.0. In the unusual case where discrete level data are not available for a particular nuclide, γ-ray emission is simulated in the same manner as for unbound nuclear states.
The TALYS structure data files are used with minor reformatting for all nuclides except 40 K. To ensure consistency with the level assignments made in section III A, I prepared an original decay scheme for 40 K using the experimental data in ref. [64], the TALYS 1.6 structure file, and (where needed) estimated γ-ray branching ratios computed using the strength function defined in section IV B 3.

B. Unbound states: statistical emission
The MARLEY approach to modeling de-excitations of unbound nuclear states rests on the assumption of compound nucleus formation: the neutrino interaction leaves the nucleus in a state of thermal equilibrium which deexcites independently of the manner in which it was formed. The number of open decay channels is taken to be large enough that competition between them can be modeled statistically. That is, transitions to individual nuclear final states can be neglected in favor of averaging over many states of approximately the same energy [79]. This last assumption is not strictly true for excitation energies slightly above the lowest fragment emission threshold. In such cases it is adopted as an approximation.
Compound nucleus modeling is a key ingredient in nuclear reaction codes designed to compute nucleon-nucleus and nucleus-nucleus cross sections, such as TALYS, EMPIRE [80], CCONE [81], and CoH 3 [82]. The treatment in MARLEY uses the Hauser-Feshbach formalism [83] common to these other codes and shares many implementation details with them.

Differential decay widths
For the present application to neutrino-induced deexcitations, the physics content of the Hauser-Feshbach statistical model (HFSM) may be conveniently summarized by the differential decay widths and which describe de-excitations of a compound nuclear state via emission of a fragment a or a γ-ray, respectively.
Here the initial (final) nucleus has excitation energy E x (E x ), total spin J (J ), and parity Π (Π ); s, , and j are the spin, orbital, and total angular momentum quantum numbers of the emitted fragment; ρ i (ρ f ) is the density of nuclear levels in the vicinity of the initial (final) state; ε is the total kinetic energy of the final particles in the rest frame of the initial nucleus; and E γ is the energy of the emitted γ-ray. For decays involving emission of a fragment with parity π a , the value of Π is fixed by conservation of parity: The possible γ-ray transitions are labeled by their multipolarity λ ≥ 1 and by whether they are electric (X = E) or magnetic (X = M) in nature. These two alternatives are distinguished based on the multipolarity and the nuclear parity: The transmission coefficients T j and T Xλ quantify how likely each decay mode is to occur. The methods used to compute them are described in sections IV B 2 and IV B 3. For practical calculations, the infinite sums that appear in eqs. (37) and (38) must be truncated. Because the value of T j (T Xλ ) falls off rapidly with increasing (λ), terms beyond = λ = 5 are neglected.
The HFSM is often communicated in terms of nuclear scattering cross sections instead of decay widths. To aid the reader in connecting the expressions given here with more standard presentations (see, e.g., refs. [84,85]), a brief derivation of eq. (37) is provided in appendix A.

Fragment transmission coefficients
The fragment transmission coefficients T j used in eq. (37) are computed by solving the radial Schrödinger equation (with relativistic kinematics as recommended in ref. [86]) is the magnitude of its 3-momentum in the rest frame of the initial nucleus, and m a (M ) denotes the mass of the emitted fragment (final nucleus). The global nuclear optical potential U developed by Koning and Delaroche [86] is used in the present calculations. A full description thereof is given in appendix C. The quantity is the fragment's kinetic energy in the rest frame of the final nucleus. The label lab is applied to this variable because it also represents the laboratory-frame kinetic energy for the time-reversed process in which the fragment is absorbed to form the compound nucleus (see appendix A). The transmission coefficient T j is related to the energy-averaged S-matrix element S j via The latter quantity may be determined by comparing the full solution u j of eq. (41) to the asymptotic form which is evaluated in terms of the proton number z (Z ) of the emitted fragment (final nucleus) and the relative speed The numerical techniques used to evaluate the fragment transmission coefficients T j are documented in ref. [

Gamma-ray transmission coefficients
The γ-ray transmission coefficients T Xλ used in eq. (38) may be written in terms of a strength function f Xλ (E γ ) such that To compute γ-ray strength functions in this work, I adopt the Standard Lorentzian model from RIPL-3 [78], which is based on early studies by Brink [88] and Axel [89]. This model assumes that γ-ray emissions of type Xλ occur via de-excitation of the corresponding giant multipole resonance, which is parameterized in terms of its centroid excitation energy E X , width Γ X , and peak cross section σ X . The strength function is evaluated according

Transitions to discrete nuclear levels
The differential decay widths given in eqs. (37) and (38) are appropriate for use at high excitation energies E x where the nuclear levels may be modeled as a continuum. When a discrete level scheme is available for the final-state nuclide, MARLEY uses the excitation energy of the last tabulated level as the lower bound for the continuum. Otherwise, a continuum level density ρ f is used all the way down to the ground state (E x = 0).
Decays to discrete levels of the final-state nucleus are considered by MARLEY in terms of the HFSM partial decay widths and The symbol δ π , which enforces parity conservation, is equal to one if eq. (39) is satisfied and zero if it is not. If J + J < 1, the width Γ γ vanishes. The expressions in eqs. (50) and (51) may be derived from eqs. (37) and (38) by treating ρ f as a delta function centered on the nuclear level of interest.

Nuclear level density
In the continuum, the final level density ρ f is computed according to the Back-shifted Fermi gas model (BFM) described in appendix B. The initial level density ρ i is evaluated according to the BFM at all excitation energies. However, since the overall factor involving ρ i cancels out in the evaluation of decay branching ratios, the specific model chosen for ρ i does not have any impact on the simulation results.

V. RESULTS
In this section, the MARLEY 1.2.0 implementation of the theoretical models described above is used to obtain pre-  Figure 4 shows MARLEY predictions of the total cross section for inclusive charged-current ν e absorption on 40 Ar. The important role played by the Coulomb corrections discussed in section II B is illustrated by the different curves in the plot. The default MARLEY approach to Coulomb effects, defined in eq. (24), involves choosing the smaller of two correction factors calculated using the Fermi function and using the modified effective momentum approximation (MEMA). In fig. 4, the solid black line gives the cross section calculated using the default approach, while the dotted cyan and dashed red lines give, respectively, the corresponding cross sections obtained when the Fermi function and MEMA are used unconditionally. Applying either correction leads to an enhancement of the total cross section over the uncorrected result, which is drawn as the blue dash-dotted line. The relationships between the different approaches to Coulomb corrections in the present calculation are qualitatively similar to those seen previously using a CRPA model [35], but there are some details that are different, e.g., the cross sections calculated using the Fermi function and the MEMA intersect at a neutrino energy between 50-60 MeV, about 10 MeV lower than in the CRPA result.  B. Exclusive cross sections Figure 5 presents the first calculation at supernova energies of total cross sections for exclusive final states in the reaction 40 Ar(ν e , e − ) 40 K * . Each exclusive channel is labeled in terms of its hadronic content, but zero or more de-excitation γ-rays are allowed even when not explicitly listed. Below a neutrino energy of about 10 MeV, only transitions to bound nuclear levels are energetically possible. These de-excite via γ-ray emission. Single neutron emission becomes appreciable around 15 MeV. Although the proton (7.58 MeV) and alpha particle (6.44 MeV) separation energies for 40 K are comparable to the neutron separation energy (7.80 MeV), the Coulomb barrier experienced by these charged particles suppresses their emission relative to neutrons. Throughout the remainder of this paper, calculations of flux-averaged cross sections will be reported for two distinct sources of low-energy electron neutrinos. The first of these is the approximate supernova neutrino energy spectrum described in ref. [92]. Here the dependence on the neutrino energy E ν is expressed in terms of the mean energy E ν and a shape parameter α. Based on an analysis of a simulated supernova, the authors of ref. [92] report values of E ν = 14.1 MeV and α = 2.67 for the time-integrated ν e spectrum, which I denote by SN T . I also consider four instantaneous spectra estimated using fig. 1 of ref. [92]. These are labeled SN 1 through SN 4 in chronological order. Table III gives the values of the spectral parameters For an antimuon decaying at rest (µDAR), the ν e energy spectrum is given by [93] φ where m µ is the muon mass and the neutrino energy E ν satisfies Experimental facilities which generate large numbers of stopped muons, such as the Spallation Neutron Source at Oak Ridge National Laboratory, provide a valuable opportunity to study tens-of-MeV neutrino interactions using a terrestrial source [94]. Table IV reports a wide variety of flux-averaged total cross sections for each of the electron neutrino spectra φ(E ν ) described above. For each entry in the table, the flux-averaged total cross section σ f for a specific final state f was obtained via the expression and the integrals in eqs. (56) and (57) are over the entire neutrino spectrum. Here σ L is the inclusive total cross section for transitions to a particular 40 K nuclear level L and R L (f ) is the branching ratio for the final state f when the de-excitation cascade begins at the level L. The sum in eq. (56) runs over all energeticallyaccessible nuclear levels. All quantities in eq. (56) are calculated analytically except for R L (f ), which is estimated using Monte Carlo simulations of de-excitations from every nuclear level listed in the MARLEY input file ve40ArCC Bhattacharya1998.react. The statistical uncertainty associated with each entry in table IV never exceeds ten percent and is typically much smaller. Figure 6 shows flux-averaged differential cross sections predicted by MARLEY for the lab-frame scattering cosine of the outgoing electron. The upper panel shows the total result (solid black) for the SN T ν e spectrum together with the separate contributions arising from Fermi (dashed blue) and Gamow-Teller (dotted red) transitions. The lower panel presents the same quantities for the µDAR ν e spectrum. Competition between the two linear components of the cross section gives rise to a total angular distribution that is nearly flat in both cases, with SN T being slightly forward and µDAR slightly backward.

C. Electron angle and energy distributions
A recent theoretical study [36] has pointed out that forbidden nuclear transitions, which are neglected in the present calculation, have an increasingly strong effect on the electron angular distribution as the neutrino energy grows beyond a few tens of MeV. Deviations from the linear behavior shown in fig. 6 signal the breakdown of the allowed approximation used by MARLEY. A future measurement of the 40 Ar(ν e , e − ) 40 K * angular differential cross section will thus provide a powerful constraint on the nuclear modeling needed to predict the relative contributions of the allowed and forbidden transitions. ably enhanced as one moves from the SN T spectrum to the µDAR spectrum.

D. Neutrino energy reconstruction
The energy of the incident neutrino is distributed among the final products of the 40 Ar(ν e , e − ) 40 K * reaction according to the relation where E e is the total energy of the outgoing electron and T γ , T ch , and T n are, respectively, the total kinetic energies of all de-excitation γ-rays, charged hadrons, and neutrons in the final state. The small recoil kinetic energy of the remnant nucleus is included in T ch . The change in binding energy E bind is given by the expression if electron binding energies are neglected. Here is the difference of ground-state atomic masses between the remnant nucleus (R) and the nuclear target (T = 40 Ar), m e is the electron mass, and m k (Q k ) is the mass (electric charge) of the kth nuclear de-excitation product. The sums in eq. (59) run over all particles emitted during nuclear de-excitations. The minimum possible change in binding energy, occurs for final states in which only γ-rays are emitted during nuclear de-excitations. In this case, the nuclide F = 40 K produced immediately after the primary neutrino interaction and the nuclide R remaining after deexcitations are identical.
Since an a priori correction for E min bind may be applied when reconstructing the neutrino energy for any 40 Ar(ν e , e − ) 40 K * event, eq. (58) may be usefully rewritten in the form Here I have defined the excess binding energy A useful property of the excess binding energy is that only a few discrete values of this variable are likely to occur at supernova energies. A future analysis of super-nova neutrino data may therefore attempt to correct for nonzero values of bind by tagging events in which a nucleon or a heavy fragment was emitted from the struck nucleus.
Beyond the binding energy contributions, the other terms in eq. (63) vary in the degree to which they may be reconstructed by a detector. In a liquid argon time projection chamber (LArTPC), the primary electron will produce a cm-scale ionization track which may be used to determine its energy and direction. De-excitation γrays will produce isolated small energy depositions within several tens of cm of the interaction vertex, primarily via Compton scattering on atomic electrons. Reconstruction of both of these features for supernova neutrino interactions is considered in ref. [32], with the conclusion that the energy associated with each can largely be recovered under realistic detector performance assumptions. Neutron tagging and calorimetry, on the other hand, were found to be far more challenging.
Low-energy charged nuclear fragments, such as protons and alpha particles, may also produce observable ionizations in a LArTPC. A key challenge for identifying the activity induced by these particles is that, at the energies relevant for supernova neutrinos, charged hadrons will produce mm-scale or smaller ionization tracks. These will likely be difficult to distinguish from the longer track produced by the primary electron. However, if events involving charged nuclear fragment emission can be successfully tagged, perhaps by identifying unusually large charge deposits near the start of the primary electron track, then at least some of the charged hadron kinetic energy may be recoverable.
To assess the relative importance of the various terms on the right-hand side of eq. (63), I define several observables, all of which may be interpreted as a reconstructed neutrino energy E reco ν under different, often quite optimistic, assumptions. The simplest reconstruction method involves adding the outgoing electron's total energy to the minimum possible change in binding energy for CC ν e absorption on 40 Ar: This estimate of the neutrino energy may be refined by adding the summed energies of all de-excitation γ-rays and further refined by adding the summed kinetic energies of all final-state charged hadrons I call the last of these variables the visible energy while recognizing that low-energy neutrons may nevertheless produce some observable signals in a LArTPC. Finally, I consider three possible strategies for implementing a binding energy correction via tagging of finalstate nuclear fragments. All three involve conditionally adding one or more terms to the expression for the visible energy above. Under the assumption that de-excitation neutrons may be successfully tagged, I define the reconstructed neutrino energy E n tag ≡ E reco vis + δ n 1n bind (72) in which the symbol δ n is defined to be unity when a MARLEY event contains at least one final-state neutron and zero otherwise. Similarly, under the assumption that charged nuclear fragment emission may be successfully identified, I define E p tag ≡ E reco vis + δ ch 1p bind (73) in which δ ch is unity when a MARLEY event contains a charged nuclear fragment in the final state and zero when it does not. In an ideal scenario in which both of these tagging techniques are reliable, a still more refined estimate of the neutrino energy may be obtained via bind . (74) Figure 8 shows the MARLEY prediction for flux-averaged differential sections with respect to each of these observables. The top (bottom) panel of the figure presents results for the SN T (µDAR) energy spectrum defined earlier. A solid black line is used to draw the differential cross section with respect to the true neutrino energy, while the other line styles represent the various methods for reconstructing it. The E reco vis result is not shown in the top panel since it is difficult to distinguish from the E reco e+γ one on the scale of the plot. For both spectra studied, the agreement between the reconstructed and true neutrino energy distributions improves most dramatically as one moves from using only the primary electron (E reco e , thin dashed green) to using both the electron and the de-excitation γ-rays (E reco e+γ , densely dotted brown) in the reconstruction. Although inclusion of information about charged hadrons is also seen to be helpful, the next most important improvement comes from the inclusion of binding energy corrections related to neutron tagging (E n tag , thick dashed violet). Due to the higher mean energy of the µDAR spectrum, nuclear fragment emission becomes more important relative to SN T , and the impact of the tagging-based binding energy corrections on neutrino energy reconstruction becomes more pronounced.
To further quantify the performance of each of these energy reconstruction methods, fig. 9 reports the fractional root mean square (RMS) resolution for each definition of the reconstructed neutrino energy where the sum runs over over the N b simulated events which fell into the neutrino energy bin b of interest. The vertical line seen around 4 MeV for the E reco e curve (thin dashed green) corresponds to the MARLEY energy threshold for CC ν e absorption. Because the thirdforbidden transition between the ground states of 40 Ar (J π = 0 + ) and 40 K (J π = 4 − ) is neglected under the allowed approximation, MARLEY predicts a finite resolution for E reco e even at threshold. The considerable improvements in energy resolution seen between E reco e and E reco e+γ (densely dotted brown) and between E reco e+γ and E n tag (thick dashed violet) further highlight the conclusions mentioned above with respect to fig. 8: while ν e energy reconstruction in CC absorption on argon will benefit from increased information about any final-state particle species, determining the deexcitation γ-ray energies and tagging neutrons are both particularly impactful.

VI. SUMMARY AND CONCLUSIONS
Due to the potential for DUNE to obtain a once-in-alifetime large-statistics measurement of supernova electron neutrinos, achieving a detailed understanding of tens-of-MeV neutrino-argon interactions is an investment that may yield a high scientific return. This paper expands our ability to model these interactions by providing a first calculation of exclusive cross sections for the 40 Ar(ν e , e − ) 40 K * reaction at supernova energies. The implementation of the models underlying this calculation in the MARLEY event generator enables studies of neutrino energy reconstruction to be carried out easily. The simple approach pursued in section V D reveals the substantial role that measuring the energies of de-excitation γ-rays and (though difficult) neutron tagging may play in optimizing supernova ν e energy resolution in a future analysis by DUNE. Further insights are available by using MARLEY in conjunction with a realistic detector simulation [16,32]. Two major approximations adopted in MARLEY 1.2.0 constitute limitations on the present study that should be revisited in future research. The first of these is the allowed approximation invoked during derivation of the inclusive differential cross section in section II A. In a more detailed calculation of this cross section, the factor e iq·x(n) that appears in the nuclear matrix element from eq. (8) is expanded in a series of multipoles [95] that depend on the spherical Bessel function j J (|q| r n ), where r n is the magnitude of x(n) and J is the multipole order. Terms representing forbidden nuclear transitions (J > 0) vanish in the |q| → 0 limit imposed by the allowed approximation, but their contribution to the cross section becomes increasingly important as the momentum transfer grows. Since the centroid energy of the multipole giant resonances grows with J roughly like 41 J/A 1/3 MeV [96], the inclusion of forbidden transitions should enhance neutrino scattering to high-lying unbound nuclear states which de-excite primarily via fragment emission. The degree to which this observation affects the present results may be studied in the future by combining a more detailed calculation of the inclusive differential cross section with the MARLEY de-excitation model.
The second major approximation used in this work, which is shared by nearly all calculations of exclusive cross sections for tens-of-MeV neutrino-nucleus scattering, is the compound nucleus assumption discussed in section IV B. Further investigation, both theoretical and experimental, will be needed to clarify the degree to which direct knock-out and pre-equilibrium processes may safely be neglected in models of low-energy neutrinonucleus reactions. A key question is how the transition between the compound nucleus picture, which is standard for low-energy neutrinos, and the intranuclear cascade picture, which is commonly used in models of accelerator neutrino interactions, should be handled as a function of neutrino energy.
While the current discussion has focused specifically on the description of nuclear de-excitations following CC ν e absorption on 40 Ar, the model presented in section IV is sufficiently general that it may be applied unaltered in a variety of other contexts. A natural next step is the use of MARLEY together with an inclusive description of inelastic neutral-current scattering on argon, a process for which de-excitations provide the only experimental observables apart from nuclear recoil. While measurements of tens-of-MeV neutrino-argon cross sections must be pursued to meet the needs of the DUNE supernova neutrino program, more immediate opportunities for confronting MARLEY with data may become available if the code is used to obtain predictions for other nuclei. Nearfuture measurements that could provide a detailed test of MARLEY include studies of CC ν e absorption on carbon by JSNS 2 [33] and neutrino-induced neutron production on lead, iron, and copper by COHERENT [97]. Measurements of exclusive cross sections and decay rates for processes that are closely related to neutrino interactions, such as electron-nucleus scattering [45] and muon capture [98,99], may also provide helpful model constraints. Finally, the capabilities of MARLEY may prove useful in simulating nuclear de-excitations induced by processes beyond the Standard Model, including nucleon decay [100,101]  In this appendix, I give a brief derivation of the expression in eq. (37) for the nuclear fragment emission differential decay width of a compound nucleus. A similar approach can be used to obtain the result in eq. (38) for γ-ray emission. The argument presented here is a modern version of one originally given by Weisskopf in ref. [79].
Consider the decay process i → a + f in which a compound nucleus i emits a fragment a to become a finalstate nucleus f . Adopt the same notation as in section IV B 1: the initial (final) nucleus has spin J (J ), parity Π (Π ), and mass M (M ). The emitted fragment has spin s, orbital (total) angular momentum (j), mass m a , and parity π a . Denote the initial nuclear excitation energy by E x , and let the final nuclear excitation energy lie on the small interval Within an arbitrary volume V and in the rest frame of the initial nucleus, the number of states n a+f that may be populated by the decay is given by where and Here p is the 3-momentum of the emitted fragment and ρ f is the spin-parity dependent nuclear level density (see appendix B 2) for the final-state nucleus.
By detailed balance, the decay width Γ a+f is related to the width Γ i of the time-reversed absorption process a + f → i via is the number of states in which the compound nucleus i may be formed. The absorption width may be written as is the particle flux and σ = π (2J + 1) |p| 2 (2s + 1)(2J + 1) is the compound nucleus formation cross section. Here E a (E f ) is the total energy of the emitted fragment (final nucleus) and is the total kinetic energy of the a + f system. A derivation of the expression in eq. (A8) is given in ref. [34] [104]. Similar derivations can also be found in, e.g., refs. [84,85]. Combining the results above and summing over the allowed values of J , which satisfies the triangle relation leads immediately to eq. (37).

Appendix B: Level density model
The nuclear level density model used in the present calculations is the RIPL-3 parameterization [78] of the Back-shifted Fermi gas Model (BFM), which is based on the work of Koning, Hilaire, and Goriely [73]. The "back shift" used by this model, which accounts for nucleon pairing effects, involves replacing the nuclear excitation energy E x by an effective value U defined by where the energy shift is adjusted to fit experimental data using the empirical parameter δ. The pairing factor χ pair is defined by even-even 0 odd-A −1 odd-odd. (B3)

Total level density
Under the Back-shifted Fermi gas model (BFM), the total density of nuclear levels near excitation energy E x is given by the expression [73] ρ is the Fermi gas level density and is a correction intended to suppress the unphysical divergence of ρ F (E x ) at low excitation energies.
Although a constant value for the level density parameter a LD is sometimes used, I adopt the energy-dependent functional form [105] recommended by RIPL-3 to correct for the damping of shell effects at high excitation energies: Here δW (Z, A) is the shell correction energy,ã(A) is the asymptotic value of a LD at high excitation energies, and γ is a damping parameter that represents how quickly a LD (E x , Z, A) approachesã(A). The values of these three parameters are given by the relations where δM exp (Z, A) is the measured nuclear mass excess [106,107] for the nuclide with proton number Z and mass number A, and δM LDM (Z, A) is the corresponding prediction for the nuclear mass excess using the liquid drop model [78, p. 3164] [72,108]. This work uses the global "BFM effective" values of the empirical parameters α, β, δ, and γ 0 obtained in ref. [73] using fits to nuclear level data:

Spin dependence
The density of nuclear levels ρ(E x , J, Π) with total spin J and parity Π near excitation energy E x may be written in the form where R(E x , J) is the nuclear spin distribution and π(E x , J, Π) is the parity distribution. Under the assumption that the individual nucleon spins are pointing in random directions, it can be shown [109] that the spin distribution R(E x , J) is given by [73] R(E x , J) = 2J + 1 2 σ 2 exp − (J + 1 2 ) 2 2 σ 2 .
The spin cutoff parameter σ determines the width of R(E x , J). To calculate this parameter, I adopt the expression recommended by RIPL-3 [78] in the absence of discrete level data: Here S n is the neutron separation energy, and

Appendix C: Optical potential
For the statistical model calculations reported here and implemented in MARLEY, the global nuclear optical potential developed by Koning and Delaroche [86] has been adopted. This phenomenological potential is based on fits to nucleon-nucleus scattering data and may be written in the form Here d s ≡ j(j + 1) − ( + 1) − s(s + 1).
is the eigenvalue of the spin-orbit operator 2 · s s s for a projectile with definite total angular momentum j, orbital angular momentum , and spin s. The Coulomb potential V C is taken to be that of a uniformly-charged sphere: In the expression above, r is the radial coordinate of the projectile, Z (z) is the proton number of the target nu-cleus (projectile), e is the elementary charge, and R C is the Coulomb radius of the nucleus.

Nucleon projectiles
The volume (V ), surface (D), and spin-orbit (SO) terms of the optical potential are functions that may be expressed as the product of an energy-dependent well depth and an energy-independent radial part: Here m π is the mass of a charged pion, and the well depths V V , W V , etc. are real-valued functions of the laboratory kinetic energy ε lab of the projectile.
The radial dependence in eqs. (C4)-(C8) is given by the Woods-Saxon [111] shape with effective radius R and diffuseness parameter a. Table V lists the values of the parameters needed to compute the radially-dependent parts of the nuclear optical potential. Each effective radius R j is related to its tabulated parameter r j via R j = r j A 1/3 j ∈ {V, D, SO, C} .
Note that all parameter values listed in table V are given fm, while the expressions given in the text assume natural units ( = c = 1).
The expressions for the well depths are most conveniently written in terms of the difference between the lab-frame kinetic energy of the projectile ε lab and the nuclear Fermi energy E N F for the projectile species N ∈ {p, n} of interest: Here the Coulomb contribution to V V is given by where the symbol δ pN is defined by δ pN ≡ 0 N = n 1 N = p .
(C18) Tables 10 and 11 from ref. [86] list the parameters needed to calculate the well depths for a nucleon projectile.

Complex projectiles
To compute the nuclear optical potential for complex projectiles (A > 1), MARLEY implements a superposition model based on a recommendation by Madland [112]. It is equivalent to the default treatment used by TALYS. Under this approach, the radial optical model parameters for a projectile with mass number A and proton (neutron) number Z (N ) are computed by weighting the corresponding parameters for individual nucleons: The Coulomb radius parameter r C remains unchanged from the nucleon case. The well depths are evaluated according to the relations In the expressions above, the superscript n ( p ) denotes the value of the corresponding quantity for an individual neutron (proton), e.g., V n SO (ε lab ) is the spin-orbit well depth for a neutron projectile.