Data based two-body current contribution to neutrino-nucleus cross section

A phenomenological model of two-body current (2p2h) contribution to neutrino cross section is introduced. Predictions of the Valencia model for 2p2h are modified using recent CC0pi measurements from T2K and MINERvA experiments. Our results suggest a significant increase of the 2p2h cross section at neutrino energies bigger than 1 GeV and also a redistribution of 2p2h events as function of energy and momentum transfer. This may have a big impact on neutrino energy reconstruction in neutrino oscillation parameters.


I. INTRODUCTION
One of the most important unknowns in modeling neutrino-nucleus cross sections [2] is the size of the contribution coming from two-body current (2p2h) mechanism [3][4][5]. It is important to have a precise estimate of the fraction of events originating from this mechanism because in detectors like SuperKamiokande they cannot be distinguished from charge current quasi-elastic (CCQE) scatterings on bound nucleons ν l + n → l − + p,ν l + p → l + + n, where l is lepton's flavor, n, p are neutron and proton, respectively. This leads to a bias in the neutrino energy reconstruction [6][7][8][9][10][11][12] and strongly affects the precision of neutrino oscillation parameters measurements. Over the last decade a lot of theoretical studies were done aiming to understand the situation [1, [13][14][15][16][17][18][19][20][21]. The most reliable ab initio computations exist only on a restricted phase space and for light nuclei. At larger neutrino energies theoretical model predictions differ significantly among themselves [22].
Experimental studies focus mainly on CC0π (called also CCQE-like) measurements with the signal defined as 'no pion in the final state' [23]. Most of the CC0π events originate from the CCQE mechanism, but there is a significant contribution from the two-body current mechanism and also from pion production with consequent absorption inside nucleus. The advantage of this type of measurements comes from simplicity of the definition of experimental signal. The data analysis does not depend on uncertain predictions for the hadrons in the final state. The available theoretical models for the 2p2h contribution give predictions for the final state lepton only and modeling final state hadrons is based on approximations [6] and nucleon final state interactions effects [24].
Recent CC0π measurements were done by T2K and MINERvA experiments. In both cases results are published in a form of flux averaged double differential * tomasz.bonus@ift.uni.wroc.pl † jan.sobczyk@uwr.edu.pl cross section in kinematic variables describing final state (anti-)muon. T2K measured cross section for neutrinos and antineutrinos on hydrocarbon [25,26] and on water [27,28]. MINERvA published measurements on hydrocarbon for antineutrinos [29] and neutrinos [30]. Altogether, there is a lot of information in the experimental data that has not yet been fully explored. The most important CCQE contribution to the CC0π signal is well understood thanks to electron scattering studies. It has been established that in the typical T2K kinematical region theoretical models used in neutrino community reproduce the QE peak region quite well [31]. For the pion production and absorption there have been many studies which put a lot of constraint on them [32,33]. The most uncertain is the 2p2h contribution and a natural question arises how much can be learnt about it from the CC0π measurements.
The goal of this paper is to answer this question and as a result to propose a new phenomenological mode of 2p2h. The computations are done using a NuWro Monte Carlo event generator [34], but our procedure is quite general and can be employed in other MC generators and be used in neutrino oscillation experimental studies.
Our study is inspired by the MINERvA experiment attempt to resolve events' kinematics completely with calorimetric-type measurement of the interacting (anti-)neutrino energy [35,36]. A study done in the context of GENIE Monte Carlo (MC) generator [37] allowed to identify a kinematical region where more strength from the 2p2h mechanism is needed, relative to predictions of the Valencia theoretical model [1,38]. Contrary to the above mentioned study our work uses information contained in the final state muon only.
Our paper is organized as follows. In section II our approach is presented and the data sets used in the numerical analysis are described. Section III outlines our main results: the new model and its performance compared to the experimental data. The Section IV contains a discussion of the results and final remarks. Appendices A and B include technical details supplementary to the Section II. A simple toy model illustrates our method of analyzing the data based on a separation of the covariance matrices into shape and normalization parts. The starting point for our investigation is the Valencia model of the 2p2h contribution described in Ref.
[1] with a restriction on the values of momentum transfer q ≤ 1.2 GeV/c [38]. It is implemented in NuWro in terms of five structure functions depending on energy and momentum transfers (ω, q). The Valencia model does not provide predictions for final state nucleons and this information is added in a factorization scheme using a model proposed in Ref. [6]. The structure functions W j define double differential cross section in final state lepton kinematical variables: In the above equations E is neutrino energy, m is charged lepton mass, G F is Fermi constant, θ c -Cabibbo angle, Q 2 ≡ q 2 − ω 2 and M is nucleon mass. A sign ± in the W 3 containing term refers to neutrino/antineutrino cases. At neutrino energies in current and planned shortand long-baseline oscillation experiments the contributions from W 4 and W 5 containing terms are strongly suppressed due to presence of charged lepton mass in a multiplicative factor.
Our considerations are based on the hypothesis that the overall double differential cross section defined by the Valencia model should be scaled by an unknown function S(ω, q): and a form of S(ω, q) will be deduced from the CC0π data. Equivalently, this may be viewed as a simultaneous rescaling of all the structure functions W j by S(ω, q): Even if the assumption introduced in Eq. 4 looks general it is in fact quite restrictive. The proposed rescaling is independent on neutrino energy and is the same for both neutrinos and antineutrinos. In Sect. IV we will explain how it can be made more general and realistic. The form of the scaling function S(ω, q) will be determined by minimization of χ 2 estimator introduced in II C.

A. Data sets
We investigate information from the T2K and MIN-ERvA CC0π measurements in the balanced way. Both are done on the same target but with different beams peaked at ∼ 600 MeV for T2K and ∼ 3.5 GeV for MIN-ERvA. In the case of MINERvA we include results from neutrinos [30] and antineutrinos [29]. In the case of T2K we include neutrino and antineutrino measurements from Ref. [26]. In all the considered measurements the results for double differential cross section is reported together with the covariance matrix V j,k . In the case of T2K data we use two separate covariance matrices for neutrino and antineutrino results in the same way in which the MINERvA data is available. We disregard the T2K neutrino/antineutrino covariance matrix in order to treat both experiments in a symmetric way.
MINERvA ν µ data contains 156 2D bins. They are distributed on 2-dimensional grid of the size 12 by 13. The binning is done respectively by longitudinal (range from 1.5 to 15 GeV/c) and transversal (range from 0 to 2.5 GeV/c) components of the outgoing muon momentum. For the MINERvAν µ data the division is done by using the same kinematic variables but binning is different (for the transversal component the range is from 0 to 1.5 GeV/c) resulting in a 10 by 6 grid i.e. 60 2D bins. In the MINERvA experiment there is a limited acceptance of muons: its angle must be lower than 20 o with respect to neutrino beam.
T2K data represent double differential cross section in (anti-)muon momentum and cosine of the lepton scattering angle. The binning is the same for neutrino and antineutrino. Altogether, there are 58 2D bins in each case. Muon momentum range is from 0 to 5 GeV/c. The full range of the cosine is employed. However, in the forward muon directions the binning is much finer. All the backward directions are contained in just one cosine bin extending from -1 to 0.2.

B. NuWro
NuWro [39] is a neutrino Monte Carlo generator developed at the Wroc law University starting from 2005. It can be used for neutrino energy range from ∼ 100 MeV to ∼ 100 GeV. For neutrino-nucleon scattering NuWro uses three interaction modes: CCQE [40] (and elastic for neutral current reactions), RES [41,42] which covers a region of invariant hadronic mass W ≤ 1.6 GeV and DIS including shallow and deep inelastic processes with W > 1. 6 GeV. In the case of neutrino-nucleus scattering two new interaction modes are: coherent pion production (COH) and two-body current (2p2h).
Simulations done in this paper were done using NuWro version 19.02. Nucleus is treated as a local Fermi gas (LFG). 2p2h events were generated with Valencia model [1,38]. Final state interactions play an important role for RES events and are modelled with Oset et al model [34,43]. NuWro predictions σ model k in each bin k is a sum of three contributions Schematically, our estimator is defined as: where k, l run over bins in double differential cross sections and V I;k,l is a covariance matrix for the experiment I, I = 1, ..., 4. It turns out that the function S(ω, q) obtained by minimizing Eq. 8 leads to a drastic and clearly nonphysical reduction of the cross section far below the measured cross section in most of the bins. We recognized this behavior as a manifestation of Peelle's Pertinent Puzzle (PPP) [44]. We checked that this effect comes from both MINERvA data sets. Various remedies were proposed to deal with this problem. We decided to follow the ideas proposed in [45,46]. The overall covariance matrix is decomposed into 'shape', 'normalization' and 'mixed' parts [46], see the details in Appendix A. Our estimator for the MINERvA data is constructed as a sum of contributions from shape and normalization uncertainties. This is similar to the treatment discussed in Ref. [46]. However, while in the MiniBooNE paper only the diagonal part of the shape covariance matrix is explored we include the complete information contained there. We performed several tests of the performance of this method and results are summarized in Appendix B. The final form of our estimator is: σ model norm,l are linearly rescaled model predictions satisfying V pseudoinv shape,k,l is Moore-Penrose pseudoinverse matrix [47] to the 'shape' component of the covariance matrix. N is defined as For the details about the estimator defined in Eq. 10 see Appendices A and B. χ 2 f inal is a function of S(ω, q) and we are looking for its minimum. In the numerical computations we approximate S(ω, q) by a 2D step function i.e. by a discrete set of values S mn where m, n refer to bins in the (ω, q) plane. m, n run values 1, . . . , 24. Continuity constraints are imposed on values of S m,n which as a result cannot be changed in a completely random way, see Sec. II E.

D. Fitter
A minimum of χ 2 f inal was found using a fitter based on a concept of genetic evolution algorithm [48]. It was chosen because of its flexibility and ability to escape from local minima.
At the beginning all the matrix entries describing parameters S mn are equal one, meaning no scaling whatsoever. At every iteration the fitter produces a set of 500 matrices, called generation. In each generation the matrices are sorted according to the values of χ 2 f inal (S k mn ). 10% of best performing matrices (the smallest χ 2 f inal ) are copied to the next generation as they are.
80% of the next generation is populated with the offspring from the previous one. In order to produce the offspring, two matrices are selected at random with a probability to select a matrix S k being: where χ 2 i is a value of χ 2 f inal of i th matrix and χ 2 max is the maximal value of χ 2 f inal in the generation. From these two parents a new matrix is built. In the first step its elements are taken from either of the parents with relative probabilities proportional to those given by Eq. 14. Continuity constrain is not yet checked at this point. In the second step about 5% (an exact number is sampled from binomial distribution) of the new matrix entries are selected at random to be modified. The modification is done with 50% of probability either by multiplication factor or by addition of a number. Multiplication factor is selected from a normal distribution centered at 1.0 and with a standard deviation 1.0. Negative values are excluded. In the case of addition a number is selected from a normal distribution centered at 0 with standard deviation 0.5. After every single modification is applied it is checked if the new number satisfies the continuity constraint defined in Eq. 15. If the constraint is not satisfied, the value is changed to the biggest/lowest allowed one.
The last 10% of the new generation consists of randomly generated matrices. They are created in the following way. We start with two empty matrices, A and B. Entries of the matrix A are filled with random values selected from a uniform distribution with minimum/ maximum being the lowest/highest values out of all entries from all the previous generations. Once matrix A is constructed, its entries are checked for the continuity constraint. If a given values satisfies the constraint, it is copied to the matrix B. If the constraint is not satisfied, a minimal/maximal allowed value according to whether the constraint is broken from above or from below, is inserted to matrix B instead. The order in which the values are checked is irrelevant as they are checked within original matrix A, which remains unchanged, and the constraint is symmetrical. At the end of this procedure the matrix B is added to the built generation.
All the percentage values, population sizes etc. were optimized during trial and error process of testing the performance of the algorithm. A lower bound value of 0.1 was imposed as a lowest possible bin value to prevent vanishing cross section from any region. Too high probability of bin modification led to instability and very slow convergence. Lower percentage of matrices copied to next generation slowed process of escaping from local minima. Higher number of random matrices does not help much, as we only need an access to explore new promising regions and the process of investigating them is time consuming.

E. Continuity constraints
To prevent obtaining rescaling matrices with large differences between neighbouring bins values we added the following constraint allowing for a control of the smoothness of the final matrix: where: • α is a user given parameter with a value from the range [0;1), • α min = 1 − α, • k, l go through the 4 closest neighbours of the i, j bin.
An impact of changing the values of α on the best fit value of χ 2 f inal is shown on Fig. 1. When we weaken the continuity constraint (α → 1) the value of χ 2 f inal at the best fit point becomes smaller. The value α = 0 corresponds to no rescaling at all. The values α = 1 ensures smooth and more physical scaling without sharp and narrow peaks in neighbouring bins. A computation for each value of α was performed in 10 5 iterations. The calculations for 10 values of α on CPU with 6 cores and 12 threads (2 fits were running simultaneously and matrix multiplication was parallelized to achieve 100% of CPU utilization) takes about 4 hours. The computations were performed 10 times and it was checked that the differences between obtained values of χ 2 f inal for each α were lower than 1%. The best results for each α were chosen as the final result.
The optimal value of α is evaluated by looking at the behavior of the function defined as It has a maximum at α ≈ 0.2 and it is the value used in all further considerations.

III. RESULTS
The final result for α = 0.2 is shown in Fig. 2. There are two regions where the scaling makes the 2p2h contribution bigger. The first one is for maximal values   For the MINERvA experiment we show separately contributions from 'shape' and 'normalization', see Eq. 10. We see that the final results seem to be determined by the MINERvA neutrino results with the largest number of bins. The contribution to χ 2 f inal from the MINERvA neutrinos was reduced by a factor of 1/3 at the expense of the MINERvA antineutrino contribution which was increased. This is a signal that our model is not general enough to accomodate both neutrino and antineutrino results at different energies. Still, the overall reduction  of the χ 2 is large which means that the new model agrees with the data much better. Another observation is that our method produces an improvement for neutrinos, regardless of their energies but is less successful for antineutrinos. It is a signal that the W 3 response function (see Eq. 3) which contributes with a different sign for neutrinos and antineutrinos should be rescaled separately.
In Fig. 4 we show contributions to the cross section from 2p2h events before and after rescaling for each experiment separately. For all the experiments we observe a significant redistribution of the strength always to the region of momentum transfer ∼ 700 MeV/c and energy transfer ∼ 500 MeV. For T2K it gives rise to a completely new picture. In the case of T2K neutrinos a region with a large cross section at momentum transfer ∼ 300 MeV/c and energy transfer ∼ 100 MeV mostly disappears and similar is the case of neutrinos with smaller values of energy end momentum transfers. The region of the strongest rescaling seen in Fig. 2 is not a very relevant one for all the experiments and leaves no visible trail in Fig. 4.
Another illustration of the performance of our model is seen in Fig. 5. A few typical histograms with experimental results and errors and also model predictions without and with rescaling calculated in this paper are shown together. We see that for MINERvA the overall size of 2p2h contribution is larger than for T2K because of bigger neutrino energy. In the case of T2K the rescaling does not introduce much change. Contrary to that in the case of MINERvA results rescaling makes the overall cross section much bigger. For neutrinos a very good data/MC agreement is obtained while for antineutrinos the rescaling seems to be too strong and MC predictions exceed the data points in some bins.
As an additional test we compared the values of χ 2 without correlations. This comparison is closest to the intuitive (sometimes misleading, though) assessment 'by eye' of data/MC agreement. In Table II we see that after rescaling the overall agreement is much better and the improvement comes mostly from neutrinos. For antineutrinos the model predictions are not changed much but also slightly improved. Apart from MINERvA neutrinos the values of χ 2 after rescaling are close to the number of degrees of freedom which means that the agreement is very good.

IV. DISCUSSION AND FINAL REMARKS
In this paper we propose a procedure to construct a phenomenological model of two-body current contribution to (anti-)neutrino cross section. A universal rescaling function to be applied to the predictions of the Valencia model [1] is found. Our result is specific to carbon target and also to a selection of models used in numerical computations in NuWro.
The results shown in Fig 5 indicate that a significant redistribution of 2p2h cross section is predicted and this translates into a change of values of reconstructed neutrino energy in experiments like T2K where in the Superkamiokande detector final state nucleons are not observed.
It is interesting that at larger neutrino energies the overall 2p2h cross sections strongly exceed those of the original Valencia model and seem to be close to the predictions from Martini et al model [3,4] and also SUSAv2 model [49]. We obtained also a strong increase of the values of response functions at the boundary of the Valencia model domain i.e. close to q = 1.2 GeV/c. This may be a signal that the definition of the boundary proposed in Ref. [38] is too restrictive and should be relaxed as it is in the SUSAv2 model. In a very recent paper of the Valencia group [50] it is argued that there is a large 3p-3h contribution neglected in the original papers. This makes the overall np-nh Valencia model cross section larger and closer to our final result.
The results presented in this paper are the first step in our program of construction of the phenomenological model of 2p2h. The final goal is very involved numerically and we decided to divide it into steps. The final step is to rescale three most important response matrices in an independent way. W 3 enters the cross section formula in Eq. 3 with different signs for neutrinos and antineutrinos and the results obtained in this paper suggest that it should be scaled in a different way than others. W 2 is multiplied by a neutrino energy dependent function that takes different values in MINERvA and T2K experiments. If we allow W 1 and W 2 to be scaled independently we get extra flexibility to adjust better to both data sets. Altogether, we think that with three independent rescalings we will obtain more reduction of χ 2 f inal from all the individual experiments.
In this paper we used the default NuWro version with LFG nucleus model. However, we think that the most prospective will be to use one of the models which are known to reproduce well the QE peak resulting from one-body mechanism, see Ref. [31]. One of them, the hole spectral function (SF) approach [51], is already implemented in NuWro. With CCQE events modeled with SF and three independent rescalings of W 1,2,3 using the approach described in this paper we should obtain a realistic model of 2p2h contribution on the carbon target. This work is in progress.

Appendix A: A decomposition of covariance matrix
In [46] one can find a procedure how to decompose an arbitrary N × N covariance matrix into a sum of 'shape', 'normalization' and 'mixed' parts: where x j are results of measurements and The matrix V shape jk is singular: a vector made of N identical numbers is an eigenvector to eigenvalue 0. This makes the use of V shape jk in the definition of the modified χ 2 difficult. We propose to introduce: where V pseudo shape is a Moore-Penrose pseudoinverse [47] of V shape . It is a generalization of the definition of inverse matrix, inverse and pseudoinverse matrices coincide for nonsingular matrices.ỹ j are normalized to satisfy jỹ j = j x j . (A4) We investigated the statistical properties of the estimator defined in Eq. A3. We used the covariance matrix of the MINERvA neutrino experiment studied in this paper. We produced several throws (y 1 , ..., y N ) generated with a multivariate distribution defined by (x 1 , ..., x N ) and V . For each one we calculated a 'normalized random throws' (ỹ 1 , ...,ỹ N ) obtained by applying a normalization factor f = j xj j yj :ỹ j = f · y j . Finally, we studied a distribution of values of χ 2 shape . It has the basic features of the standard χ 2 (N − 1) distribution. The difference is that the peak is less pronounced with more probability at both smaller and larger values of the random variable. It may be difficult to infer from χ 2 shape confidence intervals but it can be used safely as an estimator in a search for best fit values.

Appendix B: A toy model
In this appendix the performance of χ 2 introduced in Sec. II C is tested with a simple toy model. Numerical values are chosen to be similar to those used in the Ref. [45].
Suppose two measurements were done with the follow- Suppose also that a theoretical model predictions contains a parameter λ the value of which we would like to estimate based on the data. The model predictions for the two measurements are assumed to be: y(λ) = 7.2 + λ · 0.795 7.2 + λ · 0.805 .
The standard χ 2 (λ) estimator is defined as It can be checked that χ 2 (λ) has a minimum atλ = 0.94 and y(λ) = 7.95 7.96 . When we compare those values with the measurements we see that we obtained a puzzling result, a manifestation of the Peelle's Pertinent Puzzle [44]. Applying the procedure outlined in the Appendix A we obtain: A pseudoinverse of V shape is: In the χ 2 shape introduced in the Appendix A there is no information about the overall normalization of data points. A remedy is to add a N term defined as Finally we define: (y j (λ) − x j )Ṽ −1,pseudo jk (y k (λ) − x k ) + N .
It can be checked thatχ 2 f inal (λ) has a minimum at λ ≈ 1.396 and y(λ) = 8.31 8.32 which is a reasonable result.