The SuSA model for neutrino oscillation experiments: from quasielastic scattering to the resonance region

High precision studies of Beyond-Standard-Model physics through accelerator-based neutrino oscillation experiments require a very accurate description of neutrino-nucleus cross sections in a broad energy region, going from quasielastic scattering up to deep inelastic scattering. In this work we focus on the following processes: quasielastic scattering, two-particle-two-hole excitations, and the excitation of the first (Delta) and second (Roper) resonances of the nucleon. The nuclear model is fully relativistic and includes both one- and two-body currents. We compare our results with recent T2K and MicroBooNE data on carbon and argon targets, and present predictions for DUNE kinematics.


Introduction
The accurate description of neutrino-nucleus cross-sections in the GeV regime is essential for the interpretation of present and future neutrino oscillation experiments, aimed at precision measurements of the neutrino properties and the search for physics beyond the Standard Model [1]. In particular, the future HyperK [2] and DUNE [3] facilities are expected to measure the leptonic CP-violating phase δ CP , which could shed light on the origin of the matter/antimatter asymmetry in the universe. Encouraging results in this direction have recently been published by the T2K collaboration [4].
The extraction of the oscillation parameters entering the neutrino mixing matrix U from the measurements of the oscillation probabilities between different flavours crucially depends on precise knowledge of the neutrino energy, which must be inferred from the kinematics of the detected particles in the final state. Detectors are made of heavy nuclei (carbon, oxygen, argon) and a large part of the systematic error in the experimental analyses comes from the modelling of neutrino-nucleus interaction. The success of future experiments relies on the ability to reduce these nuclear uncertainties in a wide energy range, from the quasielastic (QE) region, corresponding to the elastic interaction of the neutrino with a single nucleon inside the nuclear target, up to the deep inelastic scattering (DIS) domain, where the probe interacts with the constituent quarks.
Nuclear models used in this context must satisfy some basic requirements. First of all, since typical energies belong to the GeV region, they must be relativistic, or at least contain relativistic corrections. The simplest fully relativistic nuclear model is the global Relativistic Fermi Gas (RFG), which constitutes a solid basis for more sophisticated models. The RFG framework allows for an exact relativistic treatment of both currents and nuclear states, but ignores NN correlations, aside from the statistical ones embodied in the Pauli principle. A semi-phenomenological improvement in the RFG model is represented by the Super Scaling Approximation (SuSA) model, which takes into account both initial and final state interactions as extracted from the analysis of electron scattering data at different kinematics and on different nuclei [26,27]. Another fully relativistic model, the relativistic mean field (RMF), has been shown to explain, from the microscopic perspective, the basic features of the SuSA approach, and has been used to build an updated version of the model (SuSAv2) [28], which has been applied to the study of neutrino reactions in the quasielastic region [29,30].
A second important feature required from a reliable nuclear model is consistency: the different kinematic regions and elementary processes should be described within the same theoretical framework. Consistency is easily accomplished in the RFG model, but difficult to achieve in more sophisticated models. For example, the available calculations of the 2p2h response are mostly performed in the RFG framework, and combining them with other contributions, evaluated using different, although more sophisticated, nuclear models, may lead to misleading or incorrect results.
In this paper, we will focus on the QE, 2p2h and resonance regions--the latter including the first (Delta) and second (Roper) excited states of the nucleon-within the RFG and SuSA models. The results will be compared with recent neutrino data from the T2K and MicroBooNE experiments and predictions will be shown for typical DuNE kinematics.
The paper is organized as follows: in Section 2, we summarize the formalism for charged current neutrino nucleus reactions induced by one-and two-body currents. The nuclear model described in Section 3 is used to derive the results presented in Section 4, where we compare the theoretical predictions with experimental data. Finally, in Section 5 we draw our conclusions and outline the future deveopments of this research.

Charged Current Neutrino-Nucleus Interactions
Let us consider the (ν l , l − ) charged-current (CC) cross-section for the process where a neutrino with given energy E ν and momentum k hits a nucleus A and a negative charge lepton l − is detected in the final state with energy E l , momentum k and scattering angle θ l . Here, X can be any unobserved hadronic system, containing one or more knocked out nucleons, pions and other mesons, etc. The corresponding cross-section is obtained from the contraction of the leptonic and hadronic tensors. The latter encodes the full dependence on the nuclear dynamics and is defined in the target rest frame as where J µ is the weak hadronic current, |A > is the initial nuclear ground state with energy E A and |n > are all the intermediate nuclear states, of energy E n and momentum p n , accessible through the current operator. The δ functions express energy and momentum conservation, ω = E ν − E l and q = k − k being the energy and momentum transferred from the probe to the hadronic system. The double differential cross-section can be expressed as the linear combination of five response functions [26] where being G = 1.166 × 10 −11 MeV −2 the Fermi weak constant and cos θ c = 0.975 the Cabibbo angle. The coefficients V K depend only on the lepton kinematics and are defined in Reference [26], while the response functions R K ≡ R K (| q|, ω), also defined in Reference [26], depend only on the three-momentum q and energy ω transferred to the nucleus. The indices C, L, T refer to the Coulomb, longitudinal and transverse components of the leptonic and hadronic currents with respect to q. The response functions are specific components of the hadronic tensor (2), which includes both one-and twobody terms.

One-Body Hadronic Tensor
The elastic N → N and resonance production N → N * processes are induced by one-body currents and can be treated simultaneously by introducing the inelasticity parameter [26], where m and m * are the nucleon and resonance mass, respectively, and q 2 = ω 2 − q 2 the squared four-momentum transfer. In the elastic case m * = m and ρ = 1. The single-nucleon tensor can be written in the general form [31] w µν 1b where the structure functions w i and u i depend on the specific process and are evaluated starting from the transition current.
In this work, we take the first two excited states of the nucleon into account, which dominate at the kinematics we are exploring: the spin 3/2, isospin 3/2, P 33 (1232) (∆) resonance and the spin 1/2 , isospin 1/2, P 11 (1440) (Roper) resonance. Higher resonances can easily be included in the calculation.
The structure functions relative to elastic scattering and to the N → ∆ transition are given in Reference [26] and will not be repeated here.
In the Roper resonance region, the N → P 11 (1440) weak current is [32] where are the vector and axial operators. The N → P 11 transition form factors F * i and G * i are given in Appendix A.

Two-Body Hadronic Tensor
Processes induced by two-body currents correspond to the interaction of the neutrino with a pair of correlated nucleons, leading to a 2p2h final state in which two nucleons are knocked out of the nuclear ground state. The nucleon-nucleon correlations can be modelled through the exchange of a meson and the resulting meson-exchange currents (MEC) are largely dominated by the pion. The kinematical region in which such processes occur corresponds to energy transfers between the quasielastic and ∆ resonance peaks, where the MECs are known to be essential in order to describe inclusive electron scattering data [33][34][35].
The diagrams contributing to the weak pionic MEC in the vacuum are shown in Figure 1 and are usually classified as contact (a,b), pion-in-flight (c), pion pole (d,e) and ∆-MEC (f-i). Explicit expressions for the two-body tensor w µν 2b for neutrino scattering can be found in Reference [23]. h1 h2 As a next step, one needs to embed the above one-and two-body elementary tensors into a nuclear model.

The SuSA Model
The simplest approach to a fully relativistic nuclear system is represented by the Relativistic Fermi Gas (RFG) model, in which the single-nucleon wave functions are free plane waves multiplied by Dirac spinors and the only correlations are the statistical ones induced by the Pauli principle. Each nucleus is characterized by a Fermi momentum k F , usually fitted to the width of the quasielastic peak in electron scattering data.
The one-body nuclear tensor, in both the quasielastic and the resonance regions, is given by where (E, p) and (E , p + q) are the on-shell energies and momenta of the initial and final hadrons, respectively, and w µν 1b is the elementary single-nucleon tensor defined in Equation (8). In the quasielastic case, ρ = 1 and an extra θ(| p + q| − k F ) must be inserted inside the integral (18) to account for the Pauli exclusion principle.
In this model, the response functions can be evaluated analytically and can be expressed in the general form where the U K are functions that depend on the nucleon-boson vertex and incorporate corrections due to the Fermi motion, while the "superscaling" function is a universal function-valid for all the one-body responses-depending on only one scaling variable ψ(| q|, ω; k F ). The latter is a specific combination of the transferred energy and momentum given by with ξ F = (k F /m) 2 + 1 − 1; λ = ω/(2m) and κ = | q|/(2m) are dimensionless Fermi kinetic energy, energy transfer and momentum transfer, respectively. Physically, the scaling variable ψ ρ represents, in the model, the minimal kinetic energy of the initial state nucleons participating in the reaction at a given | q| and ω in a nucleus characterized by the Fermi momentum k F . The RFG model has the advantage of being relativistic, and therefore represents a suitable starting point for more sophisticated models, but it is well-known that it gives a poor description of electron scattering data. These, unlike neutrino data, are very abundant and precise and can be used as a benchmark in neutrino scattering studies. It was first suggested in Reference [26] that the scaling behaviour of (e, e ) data can also be used as an input to obtain reliable predictions for neutrino-nucleus cross-sections. This idea is at the basis of the SuSA model, which essentially amounts to replacing the RFG superscaling function (20) with a phenomenological one, f SuSA (ψ), extracted by the analysis of electron scattering data as the ratio between the double differential cross section and an appropriate single-nucleon function [36,37]. The analysis of the longitudinal quasielastic data shows that this function is very weakly dependent on the momentum transfer q, providing that the latter is high enough (namely larger than about 400 MeV/c) to allow for the impulse approximation; this property is usually referred to as scaling of first kind. Moreover, the superscaling function is almost independent of the specific nucleus for mass numbers A ranging from 4 (helium) up to 198 (gold); this is known as scaling of second kind. Superscaling is the simultaneous occurrence of the two kinds of scaling and is well respected by electron scattering data in the QEP region. Scaling violations occur in the transverse channel due to non-impulsive contributions like 2p2h excitations.
The phenomenological superscaling function f SuSA incorporates effectively NN correlations and final state interactions and gives, by construction, a good agreement with (e, e ) data in a wide range of kinematics and mass numbers. The parametrization used in this work is f SuSA (ψ) = α where the parameters are fitted to the electron scattering quasielastic world data analyzed in Reference [37,38] for all the experimentally available kinematics and nuclear targets.
Here, we use the values α = 2.9883, β = 1.9438, γ = 0.6731 and δ = 3.8538, corresponding to the fit performed in Reference [39]. Two more parameters, the Fermi momentum k F (228 MeV/c for carbon and 241 MeV/c for argon) and the energy shift E s (20 MeV), are fitted for each nucleus to the experimental width and position of the quasielastic peak [39].
In Figure 2 the RFG and SuSA scaling functions, Equations (20) and (22), are compared with the world averaged longitudinal (e,e') data 1 . This comparison clearly shows that the RFG provides a rather poor description of electron scattering data and more realistic models must be applied to neutrino oscillation analyses. Note that although the scaling function f SuSA has been extracted from quasielastic data, in the present work, we assume it to also be valid in the resonance production region. This choice is motivated, on the one hand, by the RFG result, for which the universality of the superscaling function is exactly true, and on the other by the fact that the nuclear effects embodied in f SuSA are expected to not depend too strongly on the reaction channel. The superscaling function embodies nuclear effects, which account for both initial-and final-state physics. It is reasonable to assume that the initial state physics, essentially described by the nuclear spectral function, is independent of the reaction channel. On the other hand, the final state interactions of the produced hadrons with the nuclear medium in principle distort the scaling function in a different way in each channel. However, it was shown in References [35,40] that the use of a universal scaling function in the full spectrum provides a good description of electron-scattering data in a wide kinematical range. This makes us confident that the error associated with this approximation is not too large when the model is applied to neutrino scattering. It is also worth mentioning that an alternative approach was taken in References [26,41,42], where a scaling function to be used in the ∆ resonance region, different from the quasielastic one, was extracted from electron scattering data. This method provides a phenomenological description valid at transferred energies below the ∆ peak, while, at higher ω, it fails due to the opening of other inelastic channels. Studies on the microscopic origin of the scaling function have shown that the shape and size of f SuSA can be reproduced with good accuracy by the relativistic mean field model [43]. In particular, it was shown that the high-energy asymmetric tail displayed by f SuSA can mainly be ascribed to final state interactions and cannot be reproduced if the latter are neglected (plane wave impulse approximation) or treated inconsistently with the initial state (for instance, using an optical potential). The RMF model was also exploited to construct a new version of the superscaling model (SuSAv2) [28,35], where different scaling functions are used in each channel (longitudinal, transverse and axial, isoscalar and isovector), as predicted by the model in the quasielastic region. Although the differences between SuSA and SuSAv2 are not negligible, in this paper, we stick to the original SuSA model, which employs the same scaling function in all channels and consistently treats the quasielastic and inelastic processes. Further refinements of the model will be explored in future work.
The superscaling approach described above is based on the assumption that the neutrino interacts with a single nucleon (impulse approximation) and ignores the interaction of the probe with two correlated nucleons. These processes violate scaling of both kinds [44] and obey a different scaling law, theoretically predicted in Ref. [45] and well respected by experimental (e,e') data from different nuclei [40,46]. They are added to the model within the RFG framework.
The two-body nuclear tensor corresponding to the previously introduced MEC is evaluated in the RFG model as where an integral appears over all the 2p2h excitations of the RFG with two holes ( h 1 , h 2 ) and two particles ( p 1 , p 2 ) in the final state, and w µν 2b is the elementary two-body tensor represented in Figure 1 (see Ref. [23]).
The computation of the 2p2h responses in the RFG is time consuming due to the high dimensionality (7) of the integrals. For the purpose of the present work, where an extra integral over the neutrino flux must be performed before comparing the results to the experimental data, we make use of a parametrization of the numerical results obtained in References [29,30]. This parametrization gives a very accurate representation of the exact results in a wide kinematic range (momentum transfers up to 2 GeV/c) and provides an efficient way of obtaining completely equivalent results.
Further details on the SuSA+MEC model and on the connection between electron and neutrino scattering can be found in the recent review article [27].

Results
We now present the predictions of the model introduced in the previous Section and compare them to some recent experimental data. We consider two kinds of CC ν µ -nucleus data. The first are "0π" (or "QE-like") data, where only the outgoing muon is detected and the final state does not contain pions. These data are supposed to correspond mainly to quasielastic scattering (one nucleon knockout) and to 2p2h excitations (two nucleon knockout). Note that, in the latter, the pion exchanged between the correlated nucleons is always highly virtual. The second set of data is instead of inclusive type, in the sense that, again, only the final lepton is detected but the final state can contain any unobserved hadrons (one or more nucleons, pions, other mesons). In this case, the cross-section receives contribution not only from the QE and 2p2h processes, but also from the excitation of nucleon resonances, which subsequently decay into undetected nucleons and mesons. The non-resonant meson production can also contribute to the signal, but is supposed to be less important, as suggested by the results of References [47][48][49][50], and will therefore be ignored in this work.
We first compare our results with data published by the T2K [51,52] and Micro-BooNE [53] collaborations. Although the two experiments explore similar kinematics, the T2K off-axis neutrino flux is more focused than the broader MicroBooNE flux (see Figure 3) and this may have consequences on the relative contributions of different processes. Moreover, the nuclear targets are different: carbon for T2K and argon for MicroBooNE. This also can induce differences in nuclear effects that depend on the nuclear density.
The Fermi momenta employed in this work are k F = 228 MeV/c and 241 MeV/c for carbon and argon, respectively, and the energy shift E s = 20 MeV. These values were fitted to inclusive electron scattering data in Reference [39].
Before showing the results, a comment is in order concerning Pauli blocking effects. As mentioned previously, in the low (ω, | q|) regime, where these effects come into play, approaches based on the impulse approximation like the RFG and SuSA models should be undertaken with care. Nevertheless, since neutrino data also include this region, we include Pauli blocking in the SuSA model following the procedure originally proposed in Reference [54]. This generalizes the RFG prescription | p + q| > k F , valid only for a step-like momentum distribution, and amounts to the following replacement for the superscaling function In Figure 4, we show the SuSA model predictions for the T2K double differential (ν µ , µ − ) cross-section off 12 C with no pions in the final state as a function of the muon momentum p µ , for different bins of the scattering angle θ µ . The separate QE and MEC contributions are also shown. In all cases, the contribution of MEC (2p2h excitations) is sizable and necessary in order to explain the experimental data. The agreement with the data is rather good, except for the last angular bin and low p µ , corresponding to very small values of scattering angle. This is not surprising since, at these kinematic conditions, where small values of the energy and momentum transfer play a major role, superscaling ideas are not applicable and, in general, any model based on the impulse approximation is hardly reliable. In this region, nuclear collective effects can take place and different approaches, like the one based on RPA, are more appropriate.  It should also be mentioned that the CC0π cross-section could also receive a contribution from pion production, followed by re-absorption in the nucleus, a process not included in our calculation. This would require a microscopic description of pion production and its final state interactions, which is not available at present in our phenomenological model, where FSIs are effectively absorbed into the scaling function. According to NEUT [58] and GENIE [59] Monte Carlo generators, this contribution accounts for about 10% of the neutrino measured cross-section [60]. It should be added to the theoretical calculation, or subtracted from the data, for a detailed quantitative comparison, which is beyond the scope of this work.
Having validated the QE and MEC model versus 0π data, we now compare our results with inclusive data, which also receive a contribution from inelastic channels. As previously stated, in our approach, we include the excitation of the first two nucleon resonances, the P 33 (1232) (∆) and the P 11 (1440) (Roper).
In Figure 5, we compare the SuSA predictions with the T2K inclusive double differential (ν µ , µ − ) cross-section off 12 C, displayed versus the muon momentum p µ for different bins of the scattering angle θ µ . The analysis of the separate QE, MEC, ∆ and P 11 contributions, also shown in the figure, indicates that the ∆ resonance offers a larger contribution than the MEC and is essential to explaining the data, in particular at a small p µ , whereas the contribution of the Roper resonance is totally negligible. Some disagreement with the data at large p µ is observed for the most forward bin. This might be due to the lack of higher inelasticities in the model and will be explored in future work.   Figure 5. The SuSA inclusive double differential (ν µ , µ − ) cross-section off 12 C, averaged over the T2K flux, is displayed versus the muon momentum p µ . The separate QE, MEC, ∆ and P 11 contributions are shown. Data from Reference [52].
Similar comments hold for the MicroBooNE inclusive cross-section, shown in Figure 6. The comparison with these data is important to test the model for the argon nucleus, which will be the preferred target of future experiments. With respect to the T2K case ( Figure 5), we observe a better agreement with the experimental result at high p µ and an underestimation of the data at low p µ . The former is simply due to the larger errorbars in the experimental data, whereas the latter will likely be eliminated with the inclusion of higher inelasticities, which, for MicroBooNE, are expected to play a more important role due to the broader neutrino flux. Work along these lines is in progress. As in the case of T2K, we stress that this is a preliminary work towards a more detailed and systematic comparison model/data. For this reason, we chose not to calculate any χ 2 , but to superimpose the theoretical curves to the experimental data in order to qualitatively show the successes and deficiencies of the model. A more quantitative and complete analysis will be performed in future work. . The SuSA inclusive double differential (ν µ , µ − ) cross-section off 40 Ar, averaged over the MicroBooNE flux, is displayed, versus the muon momentum p µ . The separate QE, MEC, ∆ and P 11 contributions are shown. Data from Reference [53].
Finally, in Figures 7 and 8, we present the predictions of the SuSA model for the future DUNE experiment, characterized by a higher energy and a broader flux (see Figure 3). In this case, the contribution of the ∆ resonance becomes comparable to, or even larger than, the quasielastic one and the second resonance, P 11 plays a non-negligile, although small, role.  total Figure 8. The SuSA predictions for inclusive single differential (ν µ , µ − ) cross-section off 40 Ar, averaged over the DuNE flux, are displayed versus the muon momentum p µ . The separate QE, MEC, ∆ and P 11 contributions are shown.

Conclusions
We have presented a unified treatment of the neutrino-nucleus response from the quasielastic up to the resonance region within a semi-phenomenological nuclear model (SuSA) based on the superscaling behaviour of inclusive electron scattering data. The approach is relativistic-as required by the kinematics-and, unlike the simpler relativistic Fermi gas model or other non relativistic models, it provides a good description of electron scattering data in a wide range of kinematics, a necessary test for models used in the analysis of neutrino oscillation experiments. Moreover, the model is simple enough to be implementable in Monte Carlo generators used in the experimental analyses [61].
The SuSA model has been extensively studied in past work (see [27] and references therein), with particular focus on the quasielastic and 2p2h regions. In this work, for the first time, the approach has been extended to study the first and second resonance regions, which will be of particular interest for the future high-energy experiment DUNE. The results of the model have been successfully compared with recent T2K and MicroBooNE data and predictions have been presented for DUNE.
Finally, it is worth pointing out that the contributions of heavier resonances to the nuclear responses, as well as interference effects, should be taken into account in order to achieve a better quantitative description of the inelastic region. The present work represents a first step towards this more ambitious program.