#  Semi-inclusive two-nucleon emission in (anti) neutrino CC scattering within
the relativistic mean field framework

V.L. Martinez-Consentinoa, A.M. Cantizania, J.E. Amaroa aDepartamento de
Física Atómica, Molecular y Nuclear, and Instituto de Física Teórica y
Computacional Carlos I, Universidad de Granada, Granada 18071, Spain

###### Abstract

This paper delves into the distribution of semi-inclusive events involving the
emission of two nucleons in (anti) neutrino charged-current scattering. The
analysis is conducted within the framework of relativistic mean field theory
applied to nuclear matter. To quantify the likelihood of such semi-inclusive
events occurring, we employ a relativistic model of meson-exchange currents
that aligns with the 2p2h inclusive cross-section. The outcomes are presented
in terms of one-fold and two-fold integrated semi-inclusive cross sections. To
highlight disparities among the various emission channels, including proton-
proton, neutron-proton, and neutron-neutron, we compare them against a purely
phase-space isotropic distribution within the center of mass of the two
nucleons. These comparisons reveal significant differences in the event
distributions, shedding light on the distinctive characteristics of each
channel.

##  I Introduction

The investigation of multinucleon emission reactions induced by neutrinos and
other electroweak probes has garnered significant attention in contemporary
physics Gal11 ; Mor12 ; For12 ; Alv14 ; Mos16 ; Ath23 . This heightened
interest stems from the compelling evidence of its presence within the
quasielastic (QE) peak region, where a substantial contribution of two-nucleon
emission coexists with one-nucleon knockout events, defying disentanglement in
inclusive experiments. In other words, in this context, the precise final
hadronic state is largely unknown, except for the possible absence of pions.
This intriguing phenomenon has been underscored in the analysis of neutrino
and antineutrino scattering experiments Agu10 ; Agu13 ; Fio13 ; Abe13 .

The significance of 2p2h processes in the inclusive cross section has been
underscored and substantiated through theoretical investigations by various
research groups Mar09 ; Nie11 ; Ama11 ; Lal12 . These studies employ diverse
models that incorporate various nuclear effects, such as meson-exchange
currents (MEC) with $\Delta$-isobar excitations, final-state interactions
(FSI), short-range correlations (SRC), the random-phase approximation (RPA)
and effective interactions. These models encompass a range of approaches,
including Fermi gas models, which can be either local or global and may
incorporate relativistic corrections. Other approaches, on the other hand,
utilize shell models or quantum hadrodynamics models Cuy16 ; Cuy17 ; Ber21 .
Ab initio methods also have unveiled that MEC exert a substantial influence on
the transverse response, a phenomenon consistent with the significant presence
of 2p2h excitations. However, it’s worth noting that in these calculations,
the contribution from distinct final states cannot be readily disentangled
Lov15 ; Lov20 . In the factorization scheme, which relies on the spectral
function formalism, a similar contribution attributed to 2p2h excitations has
also been identified Roc19 . However, the inclusion of these model-dependent
ingredients has resulted in noticeable disparities among the theoretical
predictions. As a consequence, numerous research endeavors have aimed to
intercompare the results to elucidate and mitigate the systematic
uncertainties inherent in neutrino data analyses Gra13 ; Mar13b ; Mar14 ;
Ama12 .

This situation has necessitated the integration of two-nucleon ejection
mechanisms into Monte Carlo (MC) neutrino event generators Cam22 . Typically,
one commences with only the inclusive cross section information in the 2p2h
channel, as provided by theoretical groups, and lacks the corresponding semi-
inclusive cross section data. In the absence of specific knowledge regarding
the distribution of the two final particles, it has become imperative to
resort to reasonable prescriptions for implementing an algorithm that
generates events with two-nucleon final states, based on given values of
momentum and energy transfer. The conventional approach Sob12 ; GENIE ; Kat15
; Dol20  involves selecting two nucleons from the Fermi sea. By ensuring
energy-momentum conservation, the four-momentum of the final hadronic state
(comprising two nucleons) can be computed. A reasonable assumption, in the
absence of more specific information, is to consider that the distribution of
the two final particles is isotropic in the center of mass (CM) frame. In
other words, within this frame, it is assumed that the two final nucleons move
in opposite directions with equal energy and opposite momentum, and the
emission angles are randomly chosen in the CM, assuming an isotropic
distribution. Once the final momenta are determined, a boost is applied to
transition to the laboratory system, obtaining the momenta of the two ejected
nucleons in this frame. Subsequently, these momenta are further propagated
within a cascade FSI model Dyt21 ; Fra22 .

Therefore, while these event generators employ theoretical or phenomenological
methods to incorporate interactions during the propagation of the final state,
it is crucial to recognize that the isotropic 2p distribution from which they
begin is an assumption that may not be entirely accurate and could introduce
uncertainties into the simulation results. The distinction between the
isotropic distribution employed in MC generators and the more realistic
distribution associated with a 2p2h microscopic model is primarily delineated
by the value of the semi-inclusive hadronic tensor for a specific
configuration of the two final particles. This crucial ingredient is absent in
most current simulations. A recent endeavor to develop a microscopic model for
calculating the distribution of two outgoing nucleons, when compared to the
isotropic distribution model utilized in NEUT Sob20 , has revealed significant
discrepancies between the two distributions. This discrepancy between the
isotropic assumption and the actual distribution of outgoing particles
highlights a crucial area of interest in the study of neutrino-induced two-
nucleon emission processes. It underscores the need for a more detailed and
precise understanding of the semi-inclusive hadronic tensor within such
processes, particularly with regard to the configurations of the final-state
nucleons.

The relevance of studying semi-inclusive reactions and implementing them in
Monte Carlo event generators is closely tied to the need for improving the
reconstruction of incident neutrino energy based on measurements of the final
state. This reconstruction would be possible if both the final leptonic and
hadronic states were known. Efforts in this direction have included recent
measurements of semi-inclusive observables, which involve the kinematics of
both the final lepton and hadron(s), conducted by the T2K, MINERvA, and
MicroBooNE collaborations Abe18 ; Lu18 ; Cai20 ; Abr20 ; Abr20b .

Simultaneously, a series of recent investigations have focused on the
theoretical Mor14 ; Van19 ; Fra20  study of semi-inclusive reactions induced
by neutrinos, where one nucleon is detected in coincidence with the lepton
Fra21 ; Bar21 ; Fra22 ; Fra22b ; Fra23 . However, the semi-inclusive one
nucleon emission also has a contribution from 2p2h processes and, in general,
multinucleon emission. In the absence of a semi-inclusive model for the 2p2h
channel, attempts to estimate this contribution have modeled it using the
Monte Carlo implementation of 2p2h in GENIE Fra22b ; Fra23 . This, in turn,
assumes that the distribution of two nucleons is isotropic in the center of
mass, neglecting again the dependence of microscopical hadronic tensor on the
semi-inclusive variables. Therefore, this shows that there is a need for
theoretical studies of microscopic semi-inclusive versions associated with the
2p2h cross sections.

In this paper, we embark on an in-depth exploration of this issue. Our
objective is to investigate the semi-inclusive cross section in two-particle
emission reactions induced by neutrinos, employing a simple yet non-trivial
microscopic model that consistently incorporates the 2p2h hadronic tensor.
This endeavor will enable us to scrutinize the distribution of the two final
nucleons and elucidate how it deviates from the isotropic distribution in the
CM frame. To achieve this, we will utilize a relativistic model of electroweak
MEC Pac03 ; Rui16  in conjunction with the Relativistic Mean Field (RMF) model
of nuclear matter. Ros80 ; Ser86 ; Weh93 .

The MEC 2p2h responses, calculated within the framework of the relativistic
Fermi gas (RFG) model, have found extensive application in numerous analyses
and calculations of inclusive neutrino and antineutrino cross sections, in
conjunction with superscaling analysis (SuSA, SuSAv2) and the spectral
function model Meg16 ; Meg16b ; Meg19 ; Iva18 , and have been implemented
within the GENIE event generator as one of the widely utilized 2p2h
parametrization Dol20 .

In the present work, we introduce several improvements to this 2p2h MEC model.

  1. 1.

First we use the RMF framework, which accounts for the effects of the
relativistic interaction between nucleons via scalar and vector potentials,
giving rise to a relativistic effective mass and vector energy.

  2. 2.

Second, we include the complete $\Delta$ propagator, encompassing both its
real and imaginary parts. The $\Delta$ propagator assumes a crucial role
within the MEC framework when a nucleon undergoes excitation to a
$\Delta$-isobar, subsequently decaying while interacting with a second
nucleon. This process constitutes a fundamental component of the mechanisms
underlying two-nucleon emission for high energy transfer.

The 2p2h MEC model, based on the RMF framework, considered in this work has
been employed to compute the inclusive 2p2h responses and cross sections in
neutrino and antineutrino charged-current (CC) scattering Mar21 ; Mar21b ;
Mar23 . This very model has also been recently utilized to analyze the
interference between the one-body current and MEC contributions within the
1p1h response Cas23 . However, the primary focus of this work lies in
exploring the semi-inclusive response associated with the emission of two
particles, as deduced from this microscopic model. In the realm of semi-
inclusive reactions, the two emitted nucleons possess known momenta in the
proton-proton (PP), neutron-proton (NP), or neutron-neutron channels. Yet, the
state of the residual nucleus remains unknown. One of the fundamental
requirements for semi-inclusive distributions is that their integral equals
the inclusive 2p2h cross section. This condition is essential to ensure the
model’s full consistency with the inclusive framework, and it constitutes one
of the fundamental tests in this study.

The study of two-nucleon semi-inclusive reactions presents an additional
challenge due to the fact that the distributions of two nucleons expand a six-
dimensional space, making it impractical to visualize results for the full
two-nucleon distribution. A commonly employed approach is to fix four or five
of the final variables, either angles or momenta, and then examine the
distributions in the remaining variables, typically one or two dimensions.
However, this method may not yield general conclusions unless the complete map
of configurations is explored, a task that is infeasible within the confines
of a few pages. In this work, as it marks our initial exploration of these
observables, we present the results in an alternative manner. Our approach
involves fixing one or two final variables and integrating over the remaining
ones. This allows us to study partial semi-inclusive cross sections, of the
one-fold or two-fold type. These one-fold and two-fold semi-inclusive cross
sections encapsulate the global information of the full 6D distribution,
effectively averaging over the remaining variables. This approach simplifies
the comparison between pure phase-space distributions and the model
incorporating the hadronic tensor. Moreover, it allows for more
straightforward conclusions to be drawn, as any differences observed in the
one-fold or two-fold distributions imply corresponding differences in the
overall distributions. By employing this strategy, we aim to provide a
comprehensive view of the semi-inclusive reactions induced by neutrinos,
shedding light on the intricacies of two-nucleon emission processes.

In Section II, we expound upon the formalism governing semi-inclusive and
inclusive two-nucleon knockout processes initiated by neutrinos. In Section
III, we delve into the method of calculation for both the one-fold and two-
fold cross sections employed in this study. This method involves partial
summations over bins within the space of exclusive variables. Section IV is
devoted to the presentation of our results, which encompass a range of lepton
kinematics relevant to neutrino scattering reactions. Finally, in Section V,
we draw our conclusions, summarizing our findings and highlighting the key
insights gleaned from this study.

##  II Formalism

###  II.1 2p2h Cross section

In this paper we consider the reactions $(\nu_{\mu},\mu^{-}N_{1}N_{2})$ and
$(\overline{\nu}_{\mu},\mu^{+}N_{1}N_{2})$, where an incident (anti) neutrino
with energy $E_{\nu}$ interact with a complex nucleus changing into a muon
with kinetic energy $T_{\mu}$ and solid angle
$\Omega_{\mu}=(\theta_{\mu},\phi_{\mu})$, with the result of the ejection of
two nucleons, $N_{1}$, $N_{2}$, with momenta $({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2})$ in the final state plus a residual nucleus. The outgoing
nucleons can be neutrons or protons, $N_{i}=N,P$. If the two nucleons are
detected in coincidence with the $\mu^{\pm}$ lepton the probability of the
reaction is determined by the semi-inclusive cross section

$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d^{3}p^{\prime}_{1}d^{3}p^{\prime}_{2}}=\sigma_{0}(k,k^{\prime})L_{\mu\nu}W^{\mu\nu}_{\rm
N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf q},\omega)$ (1)

where $L_{\mu\nu}$ is the leptonic tensor and $W^{\mu\nu}_{\rm
N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf q},\omega)$ is the
semi-inclusive hadronic tensor. Denoting the neutrino and muon four-momenta as
$k^{\nu}=(\epsilon,{\bf k})$, and $k^{\prime}{}^{\nu}=(\epsilon^{\prime},{\bf
k}^{\prime})$, respectively, the momentum transfer is ${\bf q}={\bf k}-{\bf
k}^{\prime}$ and the energy transfer is $\omega=\epsilon-\epsilon^{\prime}$.
The four-momentum transfer $Q^{\mu}=(\omega,{\bf q}$) verifies
$Q^{2}=\omega^{2}-q^{2}<0$.

The function $\sigma_{0}$ in Eq. (1) is defined by

$\sigma_{0}(k,k^{\prime})=\frac{G^{2}\cos^{2}\theta_{c}}{4\pi^{2}}\frac{k^{\prime}}{k},$
(2)

where $G=1.166\times 10^{-11}\quad\rm MeV^{-2}$ is the Fermi constant, and
$\theta_{c}$ is the Cabibbo angle, $\cos\theta_{c}=0.975$.

The leptonic tensor $L_{\mu\nu}$ appearing in Eq. (1) is given by

$\displaystyle L_{\mu\nu}$ $\displaystyle=$ $\displaystyle
k_{\mu}k^{\prime}_{\nu}+k_{\nu}k^{\prime}_{\mu}-kk^{\prime}g_{\mu\nu}\pm
i\epsilon_{\mu\nu\alpha\beta}k^{\alpha}k^{\prime\beta},$ (3)

where $+(-)$ is for neutrino (antineutrino).

The semi-inclusive hadronic tensor $W^{\mu\nu}_{\rm N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf q},\omega)$ contains the information
about the nuclear transition matrix elements of the weak current in the above
reaction.

The inclusive cross section in the two-nucleon emission channel is recovered
by integration over the two final nucleons

$\displaystyle\left(\frac{d\sigma}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}=\sigma_{0}L_{\mu\nu}\int
d^{3}p^{\prime}_{1}d^{3}p^{\prime}_{2}W^{\mu\nu}({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2},{\bf q},\omega)=\sigma_{0}L_{\mu\nu}W^{\mu\nu}_{2p2h}({\bf
q},\omega),$ (4)

with

$W^{\mu\nu}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf
q},\omega)=\sum_{N_{1}N_{2}}W^{\mu\nu}_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2},{\bf q},\omega),$ (5)

where the isospin sum is performed over final pairs (NP , PP) in the case of
neutrino scattering, and (NN, NP) for antineutrinos. The inclusive hadronic
tensor in the two-particle, two-hole (2p2h) channel is defined

$W^{\mu\nu}_{2p2h}({\bf
q},\omega)=\sum_{N_{1}N_{2}}\left(W^{\mu\nu}_{2p2h}\right)_{N_{1}N_{2}}({\bf
q},\omega)$ (6)

where the inclusive 2p2h hadronic tensor in the $(N_{1},N_{2})$ channel is

$\left(W^{\mu\nu}_{2p2h}\right)_{N_{1}N_{2}}({\bf q},\omega)=\int
d^{3}p^{\prime}_{1}d^{3}p^{\prime}_{2}W^{\mu\nu}_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf q},\omega).$ (7)

Figure 1: Feynman diagrams of meson exchange currents considered in the
present work.

In Ref. Sim17  a MEC model was developed for the inclusive 2p2h hadronic
tensor in the RFG. With independent-particle models such as the non-
interacting Fermi gas, the emission of two particles in a nuclear transition
is only possible through a two-body current operator. In this MEC model we
consider the diagrams with exchange of a pion as shown in Fig. 1. This
approach is based on the weak pion emission model of Her07 . Here we consider
the same MEC model extended Mar21  to include the RMF within the framework of
the Walecka model Ros80 ; Ser86 ; Weh93 . In this model the inclusive 2p2h
hadronic tensor in the $N_{1},N_{2}$ channel is written

$\displaystyle(W^{\mu\nu}_{\rm 2p2h})_{N_{1}N_{2}}$
$\displaystyle=\frac{V}{(2\pi)^{9}}\int
d^{3}p^{\prime}_{1}d^{3}p^{\prime}_{2}d^{3}h_{1}d^{3}h_{2}\frac{(m^{*}_{N})^{4}}{E_{1}E_{2}E^{\prime}_{1}E^{\prime}_{2}}w^{\mu\nu}_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2})\delta({\bf
p}^{\prime}_{1}+{\bf p}^{\prime}_{2}-{\bf q}-{\bf h}_{1}-{\bf h}_{2})$ (8)
$\displaystyle\times\theta(p^{\prime}_{1}-k_{F})\theta(k_{F}-h_{1})\theta(p^{\prime}_{2}-k_{F})\theta(k_{F}-h_{2})\delta(E^{\prime}_{1}+E^{\prime}_{2}-E_{1}-E_{2}-\omega),$

where $w^{\mu\nu}_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf
h}_{1},{\bf h}_{2})$ is the elementary 2p2h hadronic tensor given below, Eq.
(11), and $V/(2\pi)^{3}=Z/(\frac{8}{3}\pi k_{F}^{3})$ for symmetric nuclear
matter with Fermi momentum $k_{F}$. In Eq. (8) the variable $m_{N}^{*}$ is the
relativistic effective mass of the nucleon, the four-momenta of the final
particles and holes are $P^{\prime}_{i}{}^{\mu}=(E^{\prime}_{i},{\bf
p}^{\prime}_{i})$, and $H_{i}{}^{\mu}=(E_{i},{\bf h}_{i})$, respectively.
Momentum conservation implies $\bf
p^{\prime}_{2}=h_{1}+h_{2}+q-p^{\prime}_{1}$.

Within the RMF model of nuclear matter Ros80 ; Ser86 ; Weh93 , the nucleons
are interacting with relativistic scalar and vector potentials, denoted by
$g_{s}\phi_{0}$, and $g_{v}V_{0}$, see Ref. Ser86 . The single particle wave
functions are plane waves with momentum ${\bf p}$, and with on-shell energy
$E=\sqrt{(m_{N}^{*})^{2}+p^{2}}$. The relativistic effective mass of the
nucleon is defined by $m_{N}^{*}=m_{N}-g_{s}\phi_{0}=M^{*}m_{N}$, where
$m_{N}$ is the free nucleon mass. Additionally the nucleon acquires a positive
energy due to the repulsion by the relativistic vector potential,
$E_{v}=g_{v}V_{0}$. Thus the total nucleon energy is $E_{RMF}=E+E_{v}$. Note
that the vector energy does not appear explicitly in Eq. (8). This is because
it cancels out by subtraction of particle and hole energies in the energy
delta function. However it can be seen that the vector energy enter in the
$\Delta$ current diagrams (f–i) of Fig. 1, through the energy of the
intermediate $\Delta$ excitation (see Ref. Mar21  for details).

The nuclear states in the RMF are Slater determinants constructed with plane
waves obtained by solving the free Dirac equation with effective mass
$m_{N}^{*}$. Note that we use the same effective mass for particles and holes.
All states with momentum $h<k_{F}$ are occupied in the ground state. The 2p2h
excitations are obtained by raising two particles above the Fermi level, with
momenta $p^{\prime}_{1}$ and $p^{\prime}_{2}>k_{F}$, leaving two holes with
momenta $h_{1}$ and $h_{2}<k_{F}$. The 2p2h hadronic tensor is generated by
the neutrino interaction with the two-body MEC operator, whose matrix elements
can be written as

$\displaystyle\langle
f|J^{\mu}(Q)|i\rangle=\frac{(2\pi)^{3}}{V^{2}}\delta({\bf p}^{\prime}_{1}+{\bf
p}^{\prime}_{2}-{\bf q}-{\bf h}_{1}-{\bf
h}_{2})\frac{(m^{*}_{N})^{2}}{\sqrt{E^{\prime}_{1}E^{\prime}_{2}E_{1}E_{2}}}j^{\mu}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2}).$ (9)

The spin-isospin two-body current function $j^{\mu}({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2})$ is the sum of diagrams of Fig. 1,
described in detail in Refs. Sim17 ; Mar21 .

The elementary 2p2h hadronic tensor $w^{\mu\nu}_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2})$, where
$N_{1},N_{2}$ are the charge states of the final nucleons, corresponds to the
transition

$|{\bf h}_{1}s_{1},{\bf h}_{2}s_{2}\rangle\rightarrow|{\bf
p}^{\prime}_{1}s^{\prime}_{1},{\bf p}^{\prime}_{2}s^{\prime}_{2}\rangle,$ (10)

where $s_{i},s^{\prime}_{j}$ are the initial and final spin components. We
denote this transition, in short, by
$|1,2\rangle\rightarrow|1^{\prime},2^{\prime}\rangle$. Then the elementary
2p2h hadronic tensor is given by

$w^{\mu\nu}_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf
h}_{1},{\bf
h}_{2})=\frac{1}{2}\sum_{s_{1}s_{2}s^{\prime}_{1}s^{\prime}_{2}}j^{\mu}(1^{\prime},2^{\prime},1,2)^{*}_{A}\,j^{\nu}(1^{\prime},2^{\prime},1,2)_{A}\,,$
(11)

where the two-body current matrix element is anti-symmetrized with respect the
pp or nn pair

$j^{\mu}(1^{\prime},2^{\prime},1,2)_{A}\equiv\left\\{\begin{array}[]{ccc}j^{\mu}(1^{\prime},2^{\prime},1,2)-j^{\mu}(1^{\prime},2^{\prime},2,1),\kern
28.45274pt&\mbox{for $\nu_{\mu}NN\rightarrow\mu NP$}&\mbox{or
$\overline{\nu}_{\mu}PP\rightarrow\mu^{+}NP$},\\\
j^{\mu}(1^{\prime},2^{\prime},1,2)-j^{\mu}(2^{\prime},1^{\prime},1,2),\kern
28.45274pt&\mbox{for $\nu_{\mu}NP\rightarrow\mu PP$}&\mbox{or
$\overline{\nu}_{\mu}NP\rightarrow\mu^{+}NN$}.\\\ \end{array}\right.$ (12)

The factor $1/2$ in Eq. (11) accounts for the anti-symmetry of the two-body
wave function with respect to exchange of momenta and spin quantum numbers to
avoid double counting.

###  II.2 Probability distribution of 2p2h events

In this work we are interested in studying the probability distribution of
semi-inclusive events corresponding to CC neutrino (antineutrino) scattering,
$E_{\nu}\rightarrow(T_{\mu},\Omega_{\mu})$, with two outgoing nucleons
$N_{1},N_{2}$ with momenta ${\bf p}^{\prime}_{1}$, ${\bf p}^{\prime}_{2}$ in
the final hadronic state. In Monte Carlo event generators the semi-inclusive
event distribution is obtained by choosing two random momenta for the holes,
${\bf h}_{1}$, ${\bf h}_{2}$, below $k_{F}$. Then the total momentum and
energy of the final nucleons is computed as ${\bf p}^{\prime}={\bf h}_{1}+{\bf
h}_{2}+{\bf q}$ and $E^{\prime}=E_{1}+E_{2}+\omega$. The individual momenta
${\bf p}^{\prime}_{1}$ and ${\bf p}^{\prime}_{2}$ for an event are generated
by assuming an isotropic distribution of the two final nucleons in the CM
frame, and then they are transformed to the Lab system. The 2p2h event ${\cal
E}=({\bf h}_{1},{\bf h}_{2},{\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$ is
allowed if the conditions $p^{\prime}_{i}>k_{F}$ are verified. All 2p2h events
${\cal E}$ of this kind that are allowed for a kinematics
$(E_{\nu},T_{\mu},\Omega_{\mu})$ contribute to the inclusive cross section
with some probability that depends on the value of the hadronic tensor for the
2p2h event, $w^{\mu\nu}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf
h}_{1},{\bf h}_{2})$.

To find such probability within our MEC model we proceed as follows. We start
with the inclusive 2p2h cross section written in the way

$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}$
$\displaystyle=$ $\displaystyle\frac{V}{(2\pi)^{9}}\int
d^{3}h_{1}d^{3}h_{2}\theta(k_{F}-h_{1})\theta(k_{F}-h_{2})\frac{(m^{*}_{N})^{4}}{E_{1}E_{2}}F_{N_{1}N_{2}}({\bf
h}_{1},{\bf h}_{2}),$ (13)

where the two-hole distribution function $F_{N_{1}N_{2}}({\bf h}_{1},{\bf
h}_{2})$ determines the contribution of each pair of holes $({\bf h}_{1},{\bf
h}_{2})$. It is given by

$\displaystyle F_{N_{1}N_{2}}({\bf h}_{1},{\bf h}_{2})=\int
d^{3}p^{\prime}_{1}d^{3}p^{\prime}_{2}\frac{1}{E^{\prime}_{1}E^{\prime}_{2}}\delta(E^{\prime}_{1}+E^{\prime}_{2}-E^{\prime})\delta({\bf
p}^{\prime}_{1}+{\bf p}^{\prime}_{2}-{\bf p}^{\prime})f_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2}),$ (14)

where we have introduced the total energy and momentum in the final state,
$E^{\prime}=E_{1}+E_{2}+\omega$ and ${\bf p}^{\prime}={\bf h}_{1}+{\bf
h}_{2}+{\bf q}$, respectively. The function $f_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2})$ is

$f_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf
h}_{2})=\sigma_{0}L_{\mu\nu}w^{\mu\nu}_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2},{\bf h}_{1},{\bf
h}_{2})\theta(p^{\prime}_{1}-k_{F})\theta(p^{\prime}_{2}-k_{F}).$ (15)

Note that $P^{\prime}{}^{\mu}=(E^{\prime},{\bf p}^{\prime})$ is the total
four-momentum of the two-particle final state $|{\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2}\rangle$ in the Lab frame. The integral over $({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2})$ in Eq. (14) can be performed by going
to the center of mass system of the final nucleons. This was done in Ref.
Rui14b  to compute the 2p2h phase space in frozen approximation (for ${\bf
h}_{i}=0$). Here we extend the method to arbitrary values of hole momenta and
including the MEC dependence in the hadronic tensor.

Doubly-primed variables refer to the CM system. The total final momentum in
the CM is zero. ${\bf p}^{\prime\prime}={\bf p}_{1}^{\prime\prime}+{\bf
p}_{2}^{\prime\prime}=0$, and the total final energy $E^{\prime\prime}$ is
determined by Lorentz invariance of the squared four-momentum,
$E^{\prime\prime}=\sqrt{E^{\prime 2}-{p^{\prime}}^{2}}$.

In the CM frame the two final nucleons are going back-to-back with the same
momentum and with the same energy

$E^{\prime\prime}_{1}=E^{\prime\prime}_{2}=\frac{E^{\prime\prime}}{2}=\frac{1}{2}\sqrt{E^{\prime
2}-{p^{\prime}}^{2}}=\frac{1}{2}\sqrt{(E_{1}+E_{2}+\omega)^{2}-({\bf
h}_{1}+{\bf h}_{2}+{\bf q})^{2}}.$ (16)

The condition $E^{\prime\prime}_{1}>m_{N}$ restricts the allowed values of
$({\bf h}_{1},{\bf h}_{2})$ for which two-nucleon emission is possible.

We perform the integration with respect CM coordinates, ${\bf
p}^{\prime\prime}_{1},{\bf p}^{\prime\prime}_{2}$, using
$d^{3}p^{\prime}_{i}/E^{\prime}_{i}=d^{3}p^{\prime\prime}_{i}/E^{\prime\prime}_{i}$.
Then we integrate over ${\bf p}^{\prime\prime}_{2}$ using the delta of
momentum. The two-hole distribution function, Eq. (14), becomes

$\displaystyle F_{N_{1}N_{2}}({\bf h}_{1},{\bf h}_{2})=\int
d^{3}p^{\prime\prime}_{1}\frac{1}{(E^{\prime\prime}_{1})^{2}}\delta(2E^{\prime\prime}_{1}-E^{\prime\prime})f_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2})\theta(E^{\prime
2}-p^{\prime 2}-4m_{N}^{*2}),$ (17)

where ${\bf p}^{\prime\prime}_{2}=-{\bf p}^{\prime\prime}_{1}$. The step
function $\theta(E^{\prime 2}-p^{\prime 2}-4m_{N}^{*2})$ ensures that the
invariant mass of the final two-particle hadronic system is larger than
$2m_{N}^{*}$ and also that the total speed of the CM is $v<1$. Using
$h_{1}^{\prime\prime}dh_{1}^{\prime\prime}=E_{1}^{\prime\prime}dE_{1}^{\prime\prime}$
to integrate over the energy, $E^{\prime\prime}_{1}$ with the help of the
Dirac delta function, we obtain $E^{\prime\prime}_{1}=E^{\prime\prime}/2$, as
expected. This fixes the value of the modulus
$p^{\prime\prime}_{1}=\sqrt{(E^{\prime\prime}_{1})^{2}-(m_{N}^{*})^{2}}$,
resulting

$\displaystyle F_{N_{1}N_{2}}({\bf h}_{1},{\bf
h}_{2})=\frac{p_{1}^{\prime\prime}}{2E_{1}^{\prime\prime}}\int
d\Omega^{\prime\prime}_{1}f_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2})\theta(E^{\prime 2}-p^{\prime
2}-4m_{N}^{*2}),$ (18)

where
$d\Omega^{\prime\prime}_{1}=d\cos\theta_{1}^{\prime\prime}d\phi_{1}^{\prime\prime}$
are the angles of ${\bf p}^{\prime\prime}_{1}$ in the CM system. Note that the
function $f_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf
h}_{1},{\bf h}_{2})$ inside the integral (18) has to be evaluated for the
momenta ${\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2}$ in the Lab system. Then
once ${\bf p}^{\prime\prime}_{1},{\bf p}^{\prime\prime}_{2}$ are defined they
are boosted back to the Lab system.

The CM moves with velocity ${\bf v}={\bf p}^{\prime}/E^{\prime}$. The
direction of the velocity is defined by the unit vector ${\bf u}={\bf v}/v$.
The boost of a CM vector $(E^{\prime\prime},{\bf p}^{\prime\prime})$ to the
Lab system can be written using $\gamma=1/\sqrt{1-v^{2}}$

$\displaystyle p^{\prime}_{u}$ $\displaystyle=$
$\displaystyle\gamma(vE^{\prime\prime}+p^{\prime\prime}_{u}),$ (19)
$\displaystyle{\bf p}^{\prime}_{\perp}$ $\displaystyle=$ $\displaystyle{\bf
p}^{\prime\prime}_{\perp},$ (20)

where $p^{\prime}_{u}={\bf p}^{\prime}\cdot{\bf u}$ is the component along
${\bf u}$, and ${\bf p}^{\prime}_{\perp}$ is the invariant component
perpendicular to ${\bf u}$. Using ${\bf p}^{\prime}=p^{\prime}_{u}{\bf u}+{\bf
p}^{\prime}_{\perp}$ we can write

${\bf p}^{\prime}=\gamma(vE^{\prime\prime}+p^{\prime\prime}_{u}){\bf u}+{\bf
p}^{\prime\prime}_{\perp}=(\gamma
vE^{\prime\prime}+(\gamma-1)p^{\prime\prime}_{u}){\bf
u}+p^{\prime\prime}_{u}{\bf u}+{\bf p}^{\prime\prime}_{\perp}=(\gamma
vE^{\prime\prime}+(\gamma-1)p^{\prime\prime}_{u}){\bf u}+{\bf
p}^{\prime\prime}.$ (21)

Therefore the vectors ${\bf p}^{\prime}_{i}$ in the Lab system inside the
integral (18) can be computed as

${\bf p}^{\prime}_{i}=(\gamma vE^{\prime\prime}_{i}+(\gamma-1){\bf
p}^{\prime\prime}_{i}\cdot{\bf u}){\bf u}+{\bf p}^{\prime\prime}_{i}.$ (22)

Inserting Eq. (18) into (13) we can write the 2p2h inclusive cross section in
the $N_{1}N_{2}$ channel as

$\displaystyle\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}$
$\displaystyle=$ $\displaystyle\frac{V}{(2\pi)^{9}}\int
d^{3}h_{1}d^{3}h_{2}\theta(k_{F}-h_{1})\theta(k_{F}-h_{2})\theta(E^{\prime
2}-p^{\prime 2}-4m_{N}^{*2})$ (24)
$\displaystyle\times\frac{(m^{*}_{N})^{4}}{E_{1}E_{2}}\frac{p_{1}^{\prime\prime}}{2E_{1}^{\prime\prime}}\int
d\Omega^{\prime\prime}_{1}f_{N_{1}N_{2}}({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2},{\bf h}_{1},{\bf h}_{2}).$

This integral over eight dimensions is computed using numerical methods. Note
that with eight coordinates
$(h_{1},\theta_{1},\phi_{1},h_{2},\theta_{2},\phi_{2},\theta_{r},\phi_{r})$,
we generate the full space of exclusive events ${\cal E}=({\bf h}_{1},{\bf
h}_{2},{\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$. For clarity here we denote
by $(\theta_{r},\phi_{r})$ the two nucleon angles relative to the CM, i.e.,

$\theta_{r}=\theta^{\prime\prime}_{1},\kern
28.45274pt\phi_{r}=\phi^{\prime\prime}_{1}.$ (25)

Thus the 2p2h inclusive cross section is finally written as

$\displaystyle\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}$
$\displaystyle=$
$\displaystyle\int_{S_{F}}d^{3}h_{1}\int_{S_{F}}d^{3}h_{2}\int
d\cos\theta_{r}d\phi_{r}G_{N_{1}N_{2}}({\bf h}_{1},{\bf
h}_{2},\theta_{r},\phi_{r}),$ (26)

where the integrals over holes is performed in the Fermi sphere, $S_{F}$, that
is, for $h_{i}<k_{F}$, and the function $G_{N_{1}N_{2}}$ is

$G_{N_{1}N_{2}}({\bf h}_{1},{\bf
h}_{2},\theta_{r},\phi_{r})=\frac{V}{(2\pi)^{9}}\frac{(m^{*}_{N})^{4}}{E_{1}E_{2}}\frac{p_{1}^{\prime\prime}}{2E_{1}^{\prime\prime}}\sigma_{0}L_{\mu\nu}w^{\mu\nu}_{N_{1}N_{2}}({\bf
p}^{\prime}_{1},{\bf p}^{\prime}_{2},{\bf h}_{1},{\bf
h}_{2})\theta(p^{\prime}_{1}-k_{F})\theta(p^{\prime}_{2}-k_{F})\theta(E^{\prime
2}-p^{\prime 2}-4m_{N}^{*2}),$ (27)

with

$E^{\prime}=E_{1}+E_{2}+\omega,\kern 28.45274ptp^{\prime}=|{\bf h}_{1}+{\bf
h}_{2}+{\bf q}|,$ (28)

and the factor

$\frac{p_{1}^{\prime\prime}}{2E_{1}^{\prime\prime}}=\frac{1}{2}\sqrt{1-\frac{4(m_{N}^{*})^{2}}{(E_{1}+E_{2}+\omega)^{2}-({\bf
h}_{1}+{\bf h}_{2}+{\bf q})^{2}}}$ (29)

comes from the Jacobian of the Lorentz transformation to the CM system. Note
that the step function $\theta(E^{\prime 2}-p^{\prime 2}-4m_{N}^{*2})$ ensures
that the value of above squared-root is real. The function
$G_{N_{1}N_{2}}({\bf h}_{1},{\bf h}_{2},\theta_{r},\phi_{r})$ determines the
probability distribution of exclusive events.

###  II.3 Semi-inclusive 2p2h events

Starting from equation (26), the inclusive cross section can be calculated
numerically. We proceed by a discretization procedure by selecting a
representative set of the exclusive events that contribute to the inclusive
cross section. With this set of exclusive events all kinds of semi-inclusive
events can be generated and computed through partial sums. We must clarify the
difference between the concepts of exclusive event and semi-inclusive event in
the context of this approach.

  * •

By exclusive event we mean a set of four vectors ${\cal E}=({\bf h}_{1},{\bf
h}_{2},{\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$ that represent a particular
excitation of a 2p2h state that is compatible with conservation of energy and
momentum, and thus contributes to the integral (26). Each exclusive event, in
turn, can be expressed as a set of two hole momenta and two relative angles in
the center of mass system of the final nucleons ${\cal E}=({\bf h}_{1},{\bf
h}_{2},\theta_{r},\phi_{r})$. Note that each exclusive event carries an
implicit, fixed value of the moment and energy transfer ${\bf q},\omega$.

  * •

A semi-inclusive event refers to a fixed value of the final momenta ${\cal
E^{\prime}}=({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$ without specifying
values for the holes $({\bf h}_{1}{\bf h}_{2})$, which are not observed.
Therefore there are many pairs of nucleons, $({\bf h}_{1}{\bf h}_{2})$, that
can contribute to a given semi-inclusive event ${\cal E^{\prime}}$, which
implies in term of probability or cross section a sum (or integral) over the
contributing pairs, $({\bf h}_{1}{\bf h}_{2})$, with the given restriction
that they belong to the semi-inclusive event ${\cal E^{\prime}}$.

  * •

A partial semi-inclusive event is defined by specifying a subset of the six
coordinates $({\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$ of the final state.
In this work we will consider i) the one-fold events that correspond to fixing
one of the values
$(p^{\prime}_{1},\theta^{\prime}_{1},\phi^{\prime}_{1},p^{\prime}_{2},\theta^{\prime}_{2},\phi^{\prime}_{2})$.
And ii) the two-fold events that correspond to fixing two of them. This will
allow us to simplify the study of the semi-inclusive cross section, which
depends on six variables, attacking first the simpler problem of the analysis
of its partial integrals as will be seen below in more detail.

  * •

Finally, a relative semi-inclusive event refers to the specification of the
two relative angles $(\theta_{r},\phi_{r})$ of the final particles, in the CM
system of the final nucleons. It must be clarified that in a relative event of
this type the total momentum of the final nucleons is not constant because
many pairs of initial nucleons $({\bf h}_{1},{\bf h}_{2})$ can contribute, and
therefore the relative semi-inclusive events are not observable. But it will
be useful to analyze the distribution of these relative events in order to
study the influence of the hadronic tensor on the semi-inclusive cross
section.

The probability distribution of partial semi-inclusive events that we will
study in the next section is determined by the following one-fold semi-
inclusive cross sections

$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{i}},\kern
28.45274pt\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta^{\prime}_{i}},\kern
28.45274pt\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\phi^{\prime}_{i}},$
(30)

for $i=1,2$, obtained by integration of the six-fold semi-inclusive cross
section over five final variables. For instance

$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{1}}$
$\displaystyle=$ $\displaystyle\int
p^{\prime}_{1}{}^{2}d\Omega^{\prime}_{1}d^{3}p^{\prime}_{2}\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d^{3}p^{\prime}_{1}d^{3}p^{\prime}_{2}}.$
(31)

The two-fold semi-inclusive cross sections are the following

$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{1}dp^{\prime}_{2}},$
$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta^{\prime}_{1}d\cos\theta^{\prime}_{2}},\kern
28.45274pt\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\phi^{\prime}_{1}d\phi^{\prime}_{2}},$
(32)
$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{i}d\cos\theta^{\prime}_{j}},$
$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{i}d\phi^{\prime}_{j}},\kern
28.45274pt\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta^{\prime}_{i}d\phi^{\prime}_{j}}.$
(33)

with $i,j=1,2$.

We discretize the integration domain of Eq. (26) to generate a finite set of
bins of exclusive events. We choose the following integration variables
$(h_{1}^{3},\cos\theta_{1},\phi_{1},h_{2}^{3},\cos\theta_{2},\phi_{2},\cos\theta_{r},\phi_{r})$.
We divide the integration interval of each variable into a finite number of
sub-intervals of respective widths $(\Delta
h^{3},\Delta\cos\theta,\Delta\phi,\Delta
h^{3},\Delta\cos\theta,\Delta\phi,\Delta\cos\theta_{r},\Delta\phi_{r})$. Note
that we use the variable $h^{3}$ instead of $h$ to perform the discretization,
that goes from $0$ to $k_{F}^{3}$. This is convenient because the property
$h^{2}dh=dh^{3}/3$ to improve the precision in the numerical integral. We use
the same integration widths for the two initial nucleons. In this way the
space of exclusive events is discretized as a finite set $\\{{\cal
E}_{i}|i=1,\ldots{\cal N}\\}$. The volume element in the dicretized space of
exclusive events is

$\Delta{\cal E}=\left(\frac{1}{3}\Delta
h^{3}\Delta\cos\theta_{h}\Delta\phi_{h}\right)^{2}\Delta\cos\theta_{r}\Delta\phi_{r}.$
(34)

The inclusive cross section can be approximated as a sum over discrete
exclusive events

$\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}\simeq\sum_{i=1}^{\cal
N}G_{N_{1}N_{2}}({\cal E}_{i})\Delta{\cal E}.$ (35)

In the limit $\Delta{\cal E}\rightarrow 0$ the exact cross section is
obtained.

According to equation (35) the probability of each exclusive event, for a
given lepton kinematics $E_{\nu},T_{\mu},\Omega_{\mu}$ is given by

$P_{N_{1}N_{2}}({\cal E})=\frac{G_{N_{1}N_{2}}({\cal E})\Delta{\cal
E}}{\displaystyle\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}},$
(36)

and it is normalized to one

$\sum_{i=1}^{\cal N}P_{N_{1}N_{2}}({\cal E}_{i})=1.$ (37)

The function $G_{N_{1}N_{2}}({\cal E})$ is given by Eq. (27) were we see that
for fixed values ${\bf h}_{1}$ and ${\bf h}_{2}$ the outgoing nucleons are
distributed according to the available phase space but they are not generated
uniformly in the center of mass of the final hadronic system, because the
function $G_{N_{1}N_{2}}({\cal E})$ depends on the value of the exclusive 2p2h
hadronic tensor, $w^{\mu\nu}({\cal E})$ for the event ${\cal E}$. This
contrasts with the procedure used in most model implementations in neutrino
interaction event generators, where an isotropic distribution for the two
outgoing nucleons is assumed. This is the case for instance of the NuWro Sob12
, NEUT Sob20 , and GENIE Dol20  event generators. Something similar happens in
GiBUU implementation of the 2p2h excitations where the exclusive hadronic
tensor does not depend on the event ${\cal E}$ Lal12 . In the results section
we will study the effect of including or not the exclusive hadronic tensor in
the distribution of the final particles and the difference between the
distributions in the different charge channels (PN, PP or NN).

##  III One- and two-fold 2p2h semi-inclusive cross sections

To calculate the one- and two-fold semi-inclusive cross sections, we begin by
generating the discrete set of exclusive events corresponding to a specific
lepton kinematics. This means a uniform set of coordinates
$(h_{1}^{3},\cos\theta_{1},\phi_{1},h_{2}^{3},\cos\theta_{2},\phi_{2},\cos\theta_{r},\phi_{r})$
for holes and relative angles of the particles, providing a discreet set of
exclusive events ${\cal E}_{i}$ that generates the inclusive cross section. In
fact, the first check we make is that the sum over events of the function
$G({\cal E})$, Eq. (35), reproduces the inclusive 2p2h cross section for the
given kinematics.

Discretization in the space of exclusive events implies in particular the
discretization of semi-inclusive events $({\bf p}^{\prime}_{1},{\bf
p}^{\prime}_{2})$. We divide the intervals of possible values of the exclusive
variables into $n$ sub-intervals or bins:

$\displaystyle p^{\prime}_{i}$ $\displaystyle:$
$\displaystyle[k_{F},(p^{\prime}_{i})_{\rm
max}]=[p^{\prime}_{i}{}^{(1)},p^{\prime}_{i}{}^{(2)},\ldots,p^{\prime}_{i}{}^{(n+1)}],$
(38) $\displaystyle\cos\theta^{\prime}_{i}$ $\displaystyle:$
$\displaystyle[-1,1]=[\cos\theta^{\prime}_{i}{}^{(1)},\cos\theta^{\prime}_{i}{}^{(2)},\ldots,\cos\theta^{\prime}_{i}{}^{(n+1)}],$
(39) $\displaystyle\phi^{\prime}_{i}$ $\displaystyle:$
$\displaystyle[0,2\pi]=[\phi^{\prime}_{i}{}^{(1)},\phi^{\prime}_{i}{}^{(2)},\ldots,\phi^{\prime}_{i}{}^{(n+1)}],$
(40)

where the maximum momentum of final nucleons is $(p^{\prime}_{i})_{\rm
max}=\sqrt{(E^{\prime}_{i})^{2}_{\rm max}-(m_{N}^{*})^{2}}$ with
$(E^{\prime}_{i})_{\rm max}=E_{F}+\omega$, corresponding to a nucleon with
Fermi energy that receives all the energy transfer.

###  III.1 One-fold semi-inclusive cross sections

The one-fold semi-inclusive cross sections can be computed for each bin as
follows.

Let $X$ be one of the semi-inclusive variables
$X=p^{\prime}_{i},\cos\theta^{\prime}_{i},\phi^{\prime}_{i}$, for $i=1,2$. For
each exclusive event ${\cal E}$ we will denote by $X({\cal E})$ the
corresponding coordinate of the event. For instance if ${\cal E}=({\bf
h}_{1},{\bf h}_{2},{\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$ is an event,
then $p^{\prime}_{1}({\cal E})=p^{\prime}_{1}$. Let us now define by
$B(X^{(k)})$ the subset of events ${\cal E}$ such as $X^{(k)}\leq X({\cal
E})<X^{(k+1)}$. That is, the event ${\cal E}$ belong to the $k$-th bin of the
variable $X$, i.e, the interval $[X^{(k)},X^{(k+1)}]$

$B(X^{(k)})\equiv\left\\{{\cal E}\kern 14.22636pt\mbox{such that}\kern
14.22636ptX^{(k)}\leq X({\cal E})<X^{(k+1)}\right\\}.$ (41)

The total probability that an event belong to the bin $k$ is obtained by
summing the probabilities of all events that fall within the bin

$P_{N_{1}N_{2}}(X^{(k)})=\sum_{{\cal E}\in B(X^{(k)})}P_{N_{1}N_{2}}({\cal
E}),$ (42)

verifying

$\sum_{k}P_{N_{1}N_{2}}(X^{(k)})=\sum_{k}\sum_{{\cal E}\in
B(X^{(k)})}P_{N_{1}N_{2}}({\cal E})=1.$ (43)

Let be $\Delta X=X^{(k+1)}-X^{(k)}$ the (constant) width of the bin of the $X$
variable. We define

$S(X^{(k)})=\frac{P_{N_{1}N_{2}}(X^{(k)})}{\Delta
X}\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}.$
(44)

Using Eq (43) we obtain

$\sum_{k}S(X^{(k)})\Delta
X=\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}.$
(45)

Therefore we identify $S(X^{(k)})$ as the corresponding one-fold semi-
inclusive cross section averaged over the bin $X^{(k)}$.

$\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dX}\right|_{X^{(k)}}=S(X^{(k)})=\frac{1}{\Delta
X}\sum_{{\cal E}\in B(X^{(k)})}G_{N_{1}N_{2}}({\cal E})\Delta{\cal E}.$ (46)

By summation over semi-inclusive bins we recover the inclusive 2p2h cross
section

$\displaystyle\sum_{k}\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dX}\right|_{X^{(k)}}\Delta
X$ $\displaystyle=$
$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}.$ (47)

By this approach we compute the six one-fold 2p2h semi-inclusive cross
sections as

$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{i}}\right|_{p^{\prime(k)}_{i}}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta p^{\prime}_{i}}\sum_{{\cal
E}\in B(p^{\prime(k)}_{i})}G_{N_{1}N_{2}}({\cal E})\Delta{\cal E},$ (48)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta^{\prime}_{i}}\right|_{\cos\theta^{\prime(k)}_{i}}$
$\displaystyle=$
$\displaystyle\frac{1}{\Delta\cos\theta^{\prime}_{i}}\sum_{{\cal E}\in
B(\cos\theta^{\prime(k)}_{i})}G_{N_{1}N_{2}}({\cal E})\Delta{\cal E},$ (49)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\phi^{\prime}_{i}}\right|_{\phi^{\prime(k)}_{i}}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta\phi^{\prime}_{i}}\sum_{{\cal
E}\in B(\phi^{\prime(k)}_{i})}G_{N_{1}N_{2}}({\cal E})\Delta{\cal E},$ (50)

where the partial sums are performed over the subset of events corresponding
to the bins of $p^{\prime}_{i},\cos\theta^{\prime}_{i}$ or
$\phi^{\prime}_{i}$:

$\displaystyle B(p^{\prime}_{i}{}^{(k)})$ $\displaystyle=$
$\displaystyle\left\\{{\cal E}\kern 14.22636pt\mbox{such that}\kern
14.22636ptp^{\prime}_{1}{}^{(k)}\leq p^{\prime}_{1}({\cal
E})<p^{\prime}_{1}{}^{(k+1)}\right\\},$ (52) $\displaystyle
B(\cos\theta^{\prime}_{i}{}^{(k)})$ $\displaystyle=$
$\displaystyle\left\\{{\cal E}\kern 14.22636pt\mbox{such that}\kern
14.22636pt\cos\theta^{\prime}_{1}{}^{(k)}\leq\cos\theta^{\prime}_{1}({\cal
E})<\cos\theta^{\prime}_{1}{}^{(k+1)}\right\\},$ (53) $\displaystyle
B(\phi^{\prime}_{i}{}^{(k)})$ $\displaystyle=$ $\displaystyle\left\\{{\cal
E}\kern 14.22636pt\mbox{such that}\kern
14.22636pt\phi^{\prime}_{1}{}^{(k)}\leq\phi^{\prime}_{1}({\cal
E})<\phi^{\prime}_{1}{}^{(k+1)}\right\\}.$ (54)

###  III.2 Two-fold semi-inclusive cross sections

The above procedure to obtain the one-fold semi-inclusive cross sections is
straightforwardly extended to the case of two-fold semi-inclusive cross
sections. If $X\neq Y$ is a pair of semi-inclusive variables
$X,Y=p^{\prime}_{i},\cos\theta^{\prime}_{i},\phi^{\prime}_{i}$, for $i=1,2$,
we construct the subset $B(X^{(k)},Y^{(l)})$ of exclusive events ${\cal E}$
such that $X({\cal E})$ is inside the $k$-th bin of the variable $X$, and
$Y({\cal E})$ is inside the $l$-th bin of the variable $Y$, i.e,

$B(X^{(k)},Y^{(l)})\equiv\left\\{{\cal E}\kern 14.22636pt|\kern
14.22636ptX^{(k)}\leq X({\cal E})<X^{(k+1)}\,\mbox{and}\,Y^{(l)}\leq Y({\cal
E})<Y^{(l+1)}\right\\}.$ (56)

The total probability that an event belong to the bins $(k,l)$ of variables
$(X,Y)$ is then

$P_{N_{1}N_{2}}(X^{(k)},Y^{(l)})=\sum_{{\cal E}\in
B(X^{(k)},Y^{(l)})}P_{N_{1}N_{2}}({\cal E}).$ (57)

Again the total probability is one

$\sum_{kl}P_{N_{1}N_{2}}(X^{(k)},Y^{(l)})=\sum_{kl}\sum_{{\cal E}\in
B(X^{(k)},Y^{(l)})}P_{N_{1}N_{2}}({\cal E})=1.$ (58)

As in the previous subsection, we define

$S(X^{(k)},Y^{(l)})=\frac{P_{N_{1}N_{2}}(X^{(k)},Y^{(l)})}{\Delta X\Delta
Y}\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}.$
(59)

Using Eq (58) we obtain

$\sum_{kl}S(X^{(k)},Y^{(l)})\Delta X\Delta
Y=\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}.$
(60)

Therefore we identify $S(X^{(k)},Y^{(l)})$ as the corresponding two-fold semi-
inclusive cross section averaged over the the bins $X^{(k)},Y^{(l)}$

$\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dXdY}\right|_{X^{(k)},Y^{(l)}}=S(X^{(k)},Y^{(l)})=\frac{1}{\Delta
X\Delta Y}\sum_{{\cal E}\in B(X^{(k)},Y^{(l)})}G_{N_{1}N_{2}}({\cal
E})\Delta{\cal E}.$ (61)

Is easy to check that the summation over two-fold semi-inclusive bins, gives
again the inclusive 2p2h cross section

$\displaystyle\sum_{kl}\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dXdY}\right|_{X^{(k)},Y^{(l)}}\Delta
X\Delta Y$ $\displaystyle=$
$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}.$ (62)

Therefore the probability distribution of the two-fold events is related to
the two-fold semi-inclusive cross section by:

$\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dXdY}\right|_{X^{(k)},Y^{(l)}}=\frac{P_{N_{1}N_{2}}(X^{(k)},Y^{(l)})}{\Delta
X\Delta
Y}\left(\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}.$
(63)

Using Eq. (61), the two-fold 2p2h semi-inclusive cross sections are computed
as

$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{1}dp^{\prime}_{2}}\right|_{p^{\prime(k)}_{i},p^{\prime(l)}_{j}}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta p^{\prime}_{1}\Delta
p^{\prime}_{2}}\sum_{{\cal E}\in
B(p^{\prime(k)}_{1},p^{\prime(l)}_{2})}G_{N_{1}N_{2}}({\cal E})\Delta{\cal
E},$ (64)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta^{\prime}_{1}d\cos\theta^{\prime}_{2}}\right|_{\cos\theta^{\prime(k)}_{1},\cos\theta^{\prime(l)}_{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\Delta\cos\theta^{\prime}_{1}\,\Delta\cos\theta^{\prime}_{2}}\sum_{{\cal
E}\in
B(\cos\theta^{\prime(k)}_{1},\cos\theta^{\prime(l)}_{2})}G_{N_{1}N_{2}}({\cal
E})\Delta{\cal E},$ (65)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\phi^{\prime}_{1}d\phi^{\prime}_{2}}\right|_{\phi^{\prime(k)}_{1},\phi^{\prime(l)}_{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\Delta\phi^{\prime}_{1}\,\Delta\phi^{\prime}_{2}}\sum_{{\cal
E}\in B(\phi^{\prime(k)}_{1},\phi^{\prime(l)}_{2})}G_{N_{1}N_{2}}({\cal
E})\Delta{\cal E},$ (66)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{i}d\cos\theta^{\prime}_{j}}\right|_{p^{\prime(k)}_{i},\cos\theta^{\prime(l)}_{j}}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta
p^{\prime}_{i}\,\Delta\cos\theta^{\prime}_{j}}\sum_{{\cal E}\in
B(p^{\prime(k)}_{i},\cos\theta^{\prime(l)}_{j})}G_{N_{1}N_{2}}({\cal
E})\Delta{\cal E},$ (67)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}dp^{\prime}_{i}d\phi^{\prime}_{j}}\right|_{p^{\prime(k)}_{i},\phi^{\prime(l)}_{j}}$
$\displaystyle=$ $\displaystyle\frac{1}{\Delta
p^{\prime}_{i}\,\Delta\phi^{\prime}_{j}}\sum_{{\cal E}\in
B(p^{\prime(k)}_{i},\phi^{\prime(l)}_{j})}G_{N_{1}N_{2}}({\cal E})\Delta{\cal
E},$ (68)
$\displaystyle\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta^{\prime}_{i}d\phi^{\prime}_{j}}\right|_{\cos\theta^{\prime(k)}_{i},\phi^{\prime(l)}_{j}}$
$\displaystyle=$
$\displaystyle\frac{1}{\Delta\cos\theta^{\prime}_{i}\,\Delta\phi^{\prime}_{j}}\sum_{{\cal
E}\in B(\cos\theta^{\prime(k)}_{i},\phi^{\prime(l)}_{j})}G_{N_{1}N_{2}}({\cal
E})\Delta{\cal E}.$ (69)

###  III.3 Relative semi-inclusive cross sections

Figure 2:  Kinematics for 2p2h semi-inclusive events used in this work. The
incident neutrino defines the $z$ axis, and the muon direction define the
scattering plane. The two final nucleon momenta define with the z-axis the two
reaction planes with angles $\phi^{\prime}_{1}$ and $\phi^{\prime}_{2}$ with
respect to the scattering plane.

Guided by the above procedure we can define the two-fold, semi-inclusive cross
section for fixed relative angles of the final nucleons. We already mentioned
that said cross section is not measurable, but mathematically it is well
defined and will allow us to shed additional light on the effect of the
hadronic tensor on the distribution of semi-inclusive events, which is the
goal of this work.

Therefore, to proceed with the definition, we will consider the variables
$\cos\theta_{r}$, $\phi_{r}$ corresponding to the relative angles of each
exclusive event, which will also be distributed in bins or sub-intervals of
$[-1,1]$, in the case of $\cos\theta_{r}$, and $[0,2\pi]$ in the case of the
$\phi_{r}$ variable. The subset of exclusive events with relative angles in
bins $(k,l)$ of the variables $(\cos\theta_{r},\phi_{r})$ is now

$B(\cos\theta_{r}^{(k)},\phi_{r}^{(l)})\equiv\left\\{{\cal E}\kern
14.22636pt|\kern 14.22636pt\cos\theta_{r}^{(k)}\leq\cos\theta_{r}({\cal
E})<\cos\theta_{r}^{(k+1)}\,\mbox{and}\,\phi_{r}^{(l)}\leq\phi_{r}({\cal
E})<\phi_{r}^{(l+1)}\right\\}.$ (71)

The total probability that an event belong to the bins $(k,l)$ of variables
$(\cos\theta_{r},\phi_{r})$ is then

$P_{N_{1}N_{2}}(\cos\theta_{r}^{(k)},\phi^{(l)})=\sum_{{\cal E}\in
B(\cos\theta_{r}^{(k)},\phi_{r}^{(l)})}P_{N_{1}N_{2}}({\cal E}),$ (72)

and the two-fold 2p2h semi-inclusive cross section averaged over the bins
$\cos\theta_{r}^{(k)},\phi_{r}^{(l)}$ is defined as in Eq. (61)

$\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta_{r}d\phi_{r}}\right|_{\cos\theta_{r}^{(k)},\phi_{r}^{(l)}}=\frac{1}{\Delta\cos\theta_{r}\Delta\phi_{r}}\sum_{{\cal
E}\in B(\cos\theta_{r}^{(k)},\phi_{r}^{(l)})}G_{N_{1}N_{2}}({\cal
E})\Delta{\cal E}.$ (73)

Once more we have the result

$\displaystyle\sum_{kl}\left.\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}d\cos\theta_{r}d\phi_{r}}\right|_{\cos\theta_{r}^{(k)},\phi_{r}^{(l)}}\Delta\cos\theta_{r}\Delta\phi_{r}$
$\displaystyle=$
$\displaystyle\frac{d\sigma_{N_{1}N_{2}}}{dT_{\mu}d\Omega_{\mu}}.$ (74)

##  IV results

The coordinate system and kinematics for the description of semi-inclusive
2p2h reaction is shown in Fig. 2. We choose the $z$-axis in the direction of
the incident neutrino. The scattering plane $(x,z)$ is defined by the final
muon and the initial neutrino. The transverse component of the muon momentum
with respect to the neutrino defines the $x$-direction. The directions of the
two final momenta, ${\bf p}^{\prime}_{i}$, and the $z$ axis define two
corresponding reaction planes that form angles $\phi^{\prime}_{1}$ and
$\phi^{\prime}_{2}$, respectively, with the scattering plane. The angles
between ${\bf p}^{\prime}_{i}$ and the $z$ axis are $\theta^{\prime}_{i}$.

Analyzing the semi-inclusive cross section, which is dependent upon the six
variables $(\mathbf{p}^{\prime}_{1},\mathbf{p}^{\prime}_{2})$, can be notably
intricate due to the necessity of selecting final events within a six-
dimensional space. Therefore, in this initial exploration of the two-nucleon
emission process, we will focus on the distributions of partial semi-inclusive
events—–both one-fold and two-fold. This entails the fixation of one or two
final hadronic variables, with integration being performed across the
remaining ones.

As shown in the preceding section, the probability of such partial events is
computed through a partial summation of exclusive event probabilities within a
discretized framework. The semi-inclusive variables $\mathbf{p}^{\prime}_{1}$
and $\mathbf{p}^{\prime}_{2}$ will similarly be discretely binned.

For the sake of convenience, we will further present the outcomes in terms of
probabilities concerning the distributions of feasible events of both one-fold
and two-fold nature, while maintaining a constant neutrino kinematic
configuration. The corresponding one-fold and two-fold cross-sections are
directly proportional to these probabilities, with the proportionality factor
being precisely the inclusive 2p2h cross-section. This value, established for
the lepton kinematics, is divided by the bin volume, as indicated in Equations
(44) and (63). In this study, we will not delve into averaging over the
neutrino flux. Instead, we will operate under the assumption of an incident
neutrino possessing a predetermined energy, yielding a muon with certain
energy and scattering angle.

Figure 3: Inclusive 2p2h neutrino cross section as a function of muon energy. Each panel corresponds to a fixed muon angle. In each panel we plot the cross section for two neutrino energies compared to the cross section averaged over the neutrino flux. The neutrino energies are those that appear in table I. The arrows indicate the muon energies in the table.  Kinematics |  $E_{\nu}$ (MeV) |  $E_{\mu}$ (MeV) | $\cos\theta_{\mu}$ |  $\theta_{\mu}$ [deg] |  $q$ [MeV/c] |  $\omega$ [MeV] |  $\theta_{q}$ [deg]
---|---|---|---|---|---|---|---
1 | 950 | 600 | 0,85 | 31,79 | 545,42 | 350 | 34,78
2 | 700 | 600 | 0,85 | 31,79 | 368,77 | 100 | 57,53
3 | 800 | 300 | 0,05 | 87,13 | 834,49 | 500 | 19,64
4 | 400 | 300 | 0,05 | 87,13 | 477,07 | 100 | 36,00
5 | 700 | 200 | $-0,75$ | 138,59 | 834,93 | 500 | 7,73
6 | 300 | 200 | $-0,75$ | 138,59 | 441,85 | 100 | 14,72
Table 1:  Selected lepton kinematics studied in this work: neutrino energy,
muon energy, scattering angle and the corresponding values of $q,\omega$ and
also the angle $\theta_{q}$ between ${\bf q}$ and the neutrino direction
($z$-axis).

We will consider the six kinematics given in table I. These kinematics have
been chosen as follows. First of all we choose three scattering angles
$\cos\theta_{\mu}=0.85$, 0.05 and $-0.75$ corresponding to three angular bins
from the MiniBoonE experiment. In Fig. 3 we plot the inclusive 2p2h cross
section averaged with the neutrino flux against the muon energy, for the three
given angles. The position of the maximum on this cross section is indicated
with the arrows in each panel, that corresponds approximately to the energy of
the muon in the kinematics of table I. In each panel of Fig. 3 we also show
the 2p2h cross section for two fixed neutrino energies. The value of the first
neutrino energy in the green curves is chosen so that the position of the
maximum of the cross section roughly coincides with that of the averaged cross
section. This is indicated by long thin arrows. These neutrino energies
corresponds to kinematics 1, 3 and 5 of tab. I,

The value of the second neutrino energy in the blue curves of fig. 3 is chosen
such that $\omega=100$ MeV at the maximum of the average cross section. This
corresponds to the kinematics 2, 4 y 6 of Tab. I. A short thick arrow shows
the contribution of this second neutrino energy to the cross section before
averaging in the flux. These neutrino energies give a smaller contribution to
the 2p2h cross section.

Once the lepton kinematics are defined, the lepton tensor can be immediately
calculated with Eq. (3) taking into account that the leptonic vectors and
momentum transfer in the coordinate system of Fig. 2 are the following

$\displaystyle k^{\mu}=(E_{\nu},{\bf k})$ $\displaystyle=$
$\displaystyle(E_{\nu},0,0,E_{\nu})$ (75) $\displaystyle k^{\nu}=(E_{\mu},{\bf
k}^{\prime})$ $\displaystyle=$
$\displaystyle(E_{\mu},k^{\prime}\sin\theta_{\mu},0,k^{\prime}\cos\theta_{\mu})$
(76) $\displaystyle{\bf q}={\bf k}-{\bf k}^{\prime}$ $\displaystyle=$
$\displaystyle(-k^{\prime}\sin\theta_{\mu},0,E_{\nu}-k^{\prime}\cos\theta_{\mu})$
(77)

As explained in the last section, for each kinematic we generate exclusive
events in terms of six coordinates
$(h_{1}^{3},\cos\theta_{1},\phi_{1},h_{2}^{3},\cos\theta_{2},\phi_{2})$ for
the two holes, and two relative angles $(\cos\theta_{r},\phi_{r})$ for the
final particles. We generate a number of $7^{6}$ hole pairs and $200^{2}$
pairs of relative angles. The total number of exclusive events, ${\cal
E}=({\bf h}_{1},{\bf h}_{2},{\bf p}^{\prime}_{1},{\bf p}^{\prime}_{2})$,
generated in this way is $N_{\rm events}=4.71\times 10^{9}$ for each one of
the kinematics of Tab. I.

Now, for each exclusive event, ${\cal E}_{i}$, we compute the exclusive
hadronic tensor $w^{\mu\nu}_{N_{1}N_{2}}({\cal E}_{i})$ for PP and NP (NN and
NP) ejection in the case of neutrino (antineutrino) scattering. After
contraction with the leptonic tensor, we build the values of the function
$G_{N_{1}N_{2}}({\cal E}_{i})$ of Eq. (27) that determine the probability
distribution of events and the semi-inclusive cross sections. We sum over the
events included in the selected bins, and compute the corresponding one-fold
and two-fold semi-inclusive cross sections averaged over the bin using the
equations (46) and (61). Dividing by the inclusive cross section, the
probability distribution of partial semi-inclusive events is obtained.

In this work we are interested in comparing with the pure phase-space
isotropic distribution of final-state nucleons in the hadronic center-of-mass
frame, similar to what is done in the Monte Carlo implementations Sob12 . This
is equivalent to neglecting the dependence of the hadronic tensor
$W^{\mu\nu}_{N_{1}N_{2}}({\cal E})$ on the exclusive event ${\cal E}$. In our
calculation we simply set

$\sigma_{0}L_{\mu\nu}W^{\mu\nu}_{N_{1}N_{2}}({\cal E})=1.$ (78)

The resulting semi-inclusive event distribution is only due to the kinematics
of 2p2h phase-space (PS) and does not depend on the current matrix elements.
In this case the probability of exclusive events is computed similarly to Eq.
(36)

$P_{\rm P.S.}({\cal E})=\frac{G_{\rm P.S.}({\cal E})}{\sum_{\cal
E^{\prime}}G_{\rm P.S.}({\cal E^{\prime}})}$ (79)

where the exclusive phase-space $G_{\rm P.S.}$ function is defined similarly
to Eq. (27)

$G_{\rm P.S.}({\bf h}_{1},{\bf
h}_{2},\theta_{r},\phi_{r})=\frac{V}{(2\pi)^{9}}\frac{(m^{*}_{N})^{4}}{E_{1}E_{2}}\frac{p_{1}^{\prime\prime}}{2E_{1}^{\prime\prime}}\theta(p^{\prime}_{1}-k_{F})\theta(p^{\prime}_{2}-k_{F})\theta(E^{\prime
2}-p^{\prime 2}-4m_{N}^{*2})$ (80)

Note that the $G_{\rm P.S.}$ function do not depend on the charge of the
nucleons, it only depends on the phase-space kinematics. Using this phase-
space distribution we can define a semi-inclusive cross section whose integral
is equal to the inclusive cross-section, similarly to Eq. (81)

$\left.\frac{d\sigma_{\rm
P.S.}}{dT_{\mu}d\Omega_{\mu}dXdY}\right|_{X^{(k)},Y^{(l)}}\equiv\frac{P_{\rm
P.S.}(X^{(k)},Y^{(l)})}{\Delta X\Delta
Y}\left(\frac{d\sigma}{dT_{\mu}d\Omega_{\mu}}\right)_{2p2h}$ (81)

Figure 4: Two-fold Semi-inclusive cross section with respect to the relative
angles in the center of mass system of the two outgoing nucleons. From top to
bottom we show the phase-space, NP and PP emission with neutrinos, and NP and
NN emission with antineutrinos. The kinematics is #1 of Tab. 1.  Figure 5: The
same as Fig. 4 for the kinematics #2 of Tab. 1.

###  IV.1 Two-fold distributions

In Figs. 4 and 5 we show our results for the two-fold semi-inclusive cross
section of the relative angles in the center of mass system of the two
outgoing nucleons, given by Eq. (73), for kinematics 1 and 2 of Table 1. In
the top panels we show the phase-space results. The NP and PP channels for
neutrino scattering, $(\nu_{\mu},\mu NP)$, and $(\nu_{\mu},\mu PP)$,
respectively, are also shown separately, as well as the NP and NN channels for
antineutrino scattering, $(\overline{\nu}_{\mu},\mu^{+}NP)$, and
$(\overline{\nu}_{\mu},\mu^{+}NN)$, respectively. All of them have been
calculated including the full 2p2h hadronic tensor (11).

To begin our analysis, we first verified that the integrals of the
distributions, as shown in Figs. 4 and 5 over the relative angles, yield the
value of the inclusive 2p2h cross section for the specified kinematics and
charge channel. This serves as a valuable test of our present calculation. In
refs. Mar21 ; Mar21b ; Mar23 ; Mar23b , inclusive 2p2h response functions were
computed within a reference frame with the momentum transfer pointing along
the $z$-axis, utilizing a seven-dimensional integration scheme. In contrast,
in the present calculation the $z$-axis is aligned with the neutrino
direction, necessitating an eight-dimensional integration. Moreover, we
compute all components of the hadronic tensor, while in Refs. Mar21 ; Mar21b ;
Mar23 ; Mar23b  only the inclusive response functions where calculated. This
cross-check reinforces the consistency of our results.

In Figs. 4 and 5 the phase-space (PS) distributions are appropriately
normalized to the total neutron-proton (NP) plus proton-proton (PP) inclusive
neutrino cross section. A close inspection reveals that indeed, in the case of
neutrino scattering, the sum of the NP and PP cross sections approximately
matches the PS cross section in the upper panel. Furthermore, it is evident
that the PP cross section is approximately seven times greater than the NP
cross section. A similar value for the PP/NP ratio was reported in Rui16b
around the $\Delta$ peak of the 2p2h inclusive response.

It is noteworthy that the NP and PP distributions differ from the phase-space
distribution. This discrepancy arises from our inclusion of the exact
dependence of the hadronic tensor on the exclusive variables of 2p2h
excitations. It is noteworthy that the PS distribution is not perfectly
uniform, primarily due to the summation over holes. Consequently, the center-
of-mass momentum is not constant, resulting in an averaged distribution over
all possible values of the center of mass contributing to the semi-inclusive
observable. According to Table I, the results of Fig. 4 correspond to neutrino
energy $E_{\nu}=950$ MeV, $q=545$ MeV/c and $\omega=350$ MeV. For this case it
is apparent that according to phase space, the most probable events cluster
around the darker continuous band show in the top panel. However, a noticeable
difference is observed compared to the neutrino PP and NP distributions, with
the differences being more pronounced in the case of NP emission.
Additionally, we note that the distribution peaks for NP and PP emission occur
at different angles. In the case of antineutrinos, distinct distributions are
also observed for NP and NN emissions, which deviate from the phase-space
predictions. These observations emphasize the importance of considering the
full dependence of the hadronic tensor on the exclusive variables.

Figure 6:  Two-fold distributions, namely $P(p_{1},p_{2})$,
$P(\cos\theta_{1},\cos\theta_{2})$, and $P(\phi_{1},\phi_{2})$, for kinematics
#1. From top to bottom we show the phase space, the neutrino NP and PP, and
the antineutrino NP and NN emission channels, respectively. The units in each
panel are such that the corresponding distribution is normalized to one.

The distributions of the relative angles also exhibit a significant dependence
on the leptonic kinematics, as can be discerned when comparing them with the
case of different kinematics, as exemplified in Fig. 5. In this scenario,
characterized by a neutrino energy of 700 MeV, momentum transfer $q=369$ MeV/c
and energy transfer $\omega=100$ MeV, notable distinctions become apparent. In
the case of neutrino scattering, the PP distribution closely aligns with the
pure phase space, a feature corroborated by the location of the darker region
where the cross section is larger. Conversely, the NP distribution deviates
considerably from both the phase space and the PP distribution, with its
maximum found at backward angles $(\theta_{r})$. In this instance, the NP
cross section is slightly smaller than that of the PP, although it remains of
the same order of magnitude.

Similar trends are observed in the case of antineutrino scattering in Fig. 5.
Turning our attention to the antineutrino NN cross section, it closely
resembles the phase space distribution, whereas the NP cross section exhibits
noticeable differences. In this case, the maximum is shifted towards forward
angles in $\theta_{r}$. These nuanced variations underscore the intricacies of
the semi-inclusive reactions and emphasize the impact of the hadronic tensor
and lepton kinematics on the semi-inclusive 2p2h cross section.

The results presented in Figures 4 and 5 are expressed in terms of absolute
cross-section values. We have observed similar trends, demonstrating
deviations of the semi-inclusive distributions from those expected for pure
phase space, for all the leptonic kinematics of Tab. I. Notably, these
differences are most pronounced for PN emission. These results provide a
partial integration of the total semi-inclusive cross section. To gain a
broader view of the entire six-dimensional landscape but in two-dimensional,
averaged sections, we will now display results for the two-fold cross sections
involving pairs of the observable semi-inclusive variables
$(p^{\prime}_{i},\cos\theta_{i},\phi_{i})$. In this context, it is less
critical to ascertain the absolute value of the cross section, as we are aware
that the normalization is uniform across all distributions, with the inclusive
cross section serving as a common reference for a given kinematic scenario.
Therefore, we will present the results in terms of probabilities for each
semi-inclusive bin, understanding that the summation over bins in these
distributions equals one. Our primary objective is to examine the differences
among various semi-inclusive charge channels and their distinctions relative
to a pure phase-space distribution.

We show a more complete example of the available two-fold distributions in
Figs. 6, 7 and 8. In the interest of brevity and due to space limitations, we
have chosen to focus on representative examples of these combinations, mindful
of not overwhelming the reader with an exhaustive display of every possible
distribution. It is essential to acknowledge that illustrating every
conceivable combination is impractical within the confines of this
presentation. However, the selected distributions provide a meaningful and
insightful glimpse into the complex multi-dimensional landscape of semi-
inclusive reactions initiated by neutrinos and antineutrinos. Together, these
figures provide a comprehensive perspective on the distributions within the
semi-inclusive charge channels, enabling a deeper understanding of the
deviations from pure phase-space models.

As a first example, in Figure 6, we present the two-fold distributions, namely
$P(p^{\prime}_{1},p^{\prime}_{2})$,
$P(\cos\theta^{\prime}_{1},\cos\theta^{\prime}_{2})$, and
$P(\phi^{\prime}_{1},\phi^{\prime}_{2})$, for both neutrino and antineutrino
scattering. These distributions provide insights into the joint probability
distributions of pairs of momenta or emission angles of the two particles
within the semi-inclusive reactions. This example vividly demonstrates the
correlations between pairs of variables and the differences introduced by
considering the semi-inclusive hadronic tensor in contrast to pure phase-space
distributions. Notably, these distributions show the asymmetries present in
the final-state momenta of protons and neutrons within the NP channel,
underscoring the influence of the hadronic tensor.

Specifically, the $P(p^{\prime}_{1},p^{\prime}_{2})$ distribution (first
column of Fig. 6) reveals a strong correlation between these two variables.
The distribution occupies a narrow band around an curve centered at
$(k_{F},k_{F})$, corresponding to the minimum values the outgoing particles
can attain. The global shape of the distribution is primarily dictated by
kinematics, an does not depend on the hadronic tensor. We observe that the
distribution is centered and exhibits a maximum around emission momenta of
$p^{\prime}_{1}=p^{\prime}_{2}\sim 550$ MeV/c for PP, NN, and the phase-space
(PS) distributions. However, the NP distribution reveals a distinct pattern.
For neutrino scattering, its maximum is centered approximately at
$p^{\prime}_{1}\sim 650$ MeV/c and $p^{\prime}_{2}=450$ MeV/c, while for
antineutrino scattering, the situation is reversed, bearing in mind that
particle 1 is a neutron and particle 2 is a proton. This clear asymmetry
underscores the dissimilarity in final proton and neutron momenta within the
NP channel, a consequence of the influence of the hadronic tensor.

In the second column of Figure 6, we turn our attention to the correlations
between the emission angles $\theta^{\prime}_{1}$ and $\theta^{\prime}_{2})$.
Here we also observe a distinctive pattern: the emission angles tend to
cluster within a relatively narrow band, forming a curve that traverses from
backward-forward $(-1,1)$ to forward-backward $(1,-1)$ emissions. Within this
band, all pairs of angles exhibit fairly similar probabilities, except in the
extremes. Additionally, asymmetries in the distributions of NP emission
concerning particles 1 and 2 are observed, even though they might not be as
evident in the graph. Meanwhile, the distributions for PP, NN, and PS exhibit
clear symmetry under the $1\leftrightarrow 2$ exchange.

In the third column of Figure 6, we finally explore the correlations between
the azimuthal emission angles $\phi^{\prime}_{1}$ and $\phi^{\prime}_{2}$,
which correspond to the angles between the emission plane and the scattering
plane (as illustrated in Fig. 2). Once again, we observe a distinct
correlation shape that traces two precise trajectories in the
$(\phi_{1},\phi_{2})$ plane, ranging from $(0,\pi)$ to $(\pi,2\pi)$ and from
$(\pi,0)$ to $(2\pi,\pi)$. It’s worth noting that in spherical coordinates,
the values 0 and $2\pi$ are identified as the same point, so these two
trajectories actually form a single path. In this case, we observe a clear
$1\leftrightarrow 2$ symmetry in the distributions for PP, NN, and the phase
space around the angles of $(\pi/2,3\pi/2)$ and $(3\pi/2,\pi/2)$. However, in
the case of NP emission, this symmetry is disrupted again due to the influence
of the hadronic tensor.

Figure 7:  Two-fold distributions, $P(p_{1},\cos\theta_{1})$,
$P(p_{1},\phi_{1})$, $P(p_{1},\cos\theta_{2})$, and $P(p_{1},\phi_{2})$,
specifically for neutrino scattering and kinematics 1. In each panel from left
to right we show the phase space, NP and PP distributions.

As a second example, depicted in Figure 7, we present the two-fold
distributions $P(p^{\prime}_{1},\cos\theta^{\prime}_{1})$,
$P(p^{\prime}_{1},\phi^{\prime}_{1})$,
$P(p^{\prime}_{1},\cos\theta^{\prime}_{2})$, and
$P(p^{\prime}_{1},\phi^{\prime}_{2})$, specifically for neutrino scattering
and kinematics 1. These distributions shed light on the correlations between
momentum of the first particle and one of the four angular variables. When
comparing the phase-space, NP, and PP distributions, it is observed that they
share a similar shape. This similarity is attributed to the fact that, as in
the previous cases, the overall features are largely determined by kinematics,
specifically energy and momentum conservation. The influence of the hadronic
tensor alters the finer details, and clear differences can be observed between
NP and PP emissions when compared to the phase-space distribution.

Figure 8: Two-fold semi-inclusive distributions,
$P(\cos\theta^{\prime}_{1},\phi^{\prime}_{1})$ and
$P(\cos\theta^{\prime}_{1},\phi^{\prime}_{2})$, for neutrino scattering and
for kinematics 1. From left to right we show the phase-space, NP and PP
emission channels, respectively.

In a final example, in Figure 8, we showcase the two-fold distributions,
$P(\cos\theta^{\prime}_{1},\phi^{\prime}_{1})$ and
$P(\cos\theta^{\prime}_{1},\phi^{\prime}_{2})$, for neutrino scattering and
kinematics 1, highlighting the joint probabilities associated with the
emission angle of particle one and one of the azimuthal angles in the semi-
inclusive reaction. Once more, we observe that the overall shapes of these
distributions are similar among the three cases: phase-space, NP, and PP
emissions. However, notable differences become evident in the detailed
behavior of $P(\cos\theta^{\prime}_{1},\phi^{\prime}_{1})$, indicating that
this distribution is particularly sensitive to the influence of the hadronic
tensor. The distribution of $\cos\theta^{\prime}_{1}$ and $\phi^{\prime}_{2}$
exhibits a pronounced peak around forward emission, and it becomes smoother
away from that region. While there are still differences between the various
emission channels, these differences are not very noticeable in the figure.

Figure 9:  Two-fold distribution $P(p^{\prime}_{1},p^{\prime}_{2})$ for semi-
inclusive 2p2h neutrino scattering. From top to bottom, the rows refer to the
kinematics 2–6 of table I, Form left to right, the columns refer to the phase-
space, NP and PP distribution, respectively.

To conclude the discussion of two-fold distributions in semi-inclusive two-
nucleon emission, Figure 9 provides a comprehensive comparison of the
distribution $P(p^{\prime}_{1},p^{\prime}_{2})$. This comparison covers
kinematics 2 to 6 (rows in Fig. 9) for neutrino scattering, as outlined in
Table I. The three columns, from left to right, correspond to phase-space, NP,
and PP distributions, respectively. In this analysis, a clear correlation
between the two momenta remains evident across all cases. In the instances of
phase-space and PP emissions, the distributions exhibit clear symmetry with
respect to $1\leftrightarrow$ exchange. Then $p^{\prime}_{1}$ equals
$p^{\prime}_{2}$. at the maximum of the distribution. Conversely, the NP
distribution displays a clear asymmetry under $1\leftrightarrow$ exchange, and
$p^{\prime}_{1}\neq p^{\prime}_{2}$ at the maximum. These result suggest that
it is more likely for the neutron to possess more energy than the proton. The
consistent trend of distinct distributions of momentum pairs between the NP
and PP channels is observed across all the kinematic configurations analyzed
in this study. This consistency underscores the necessity of including the
hadronic tensor to accurately account for these differences.

Figure 10:  One-fold distributions for semi-inclusive two-nucleon emission
induced by neutrinos for kinematics 1\. In each panel, we compare the
distributions for phase-space (PS), NP, and PP emissions.

###  IV.2 One-fold distributions

Now we discuss the results for the one-fold distributions in semi-inclusive
two-nucleon emission induced by neutrinos (antineutrinos). Note that these
one-fold cross sections also contribute to semi-inclusive one-nucleon emission
when combined with the results obtained from a model of one-nucleon knockout.

In Figure 10, we present the distributions specifically for neutrino
scattering, focusing on kinematics 1 as outlined in Table I. In each panel, we
display the one-fold distributions of the variables
$p^{\prime}_{1},p^{\prime}_{2},\cos\theta^{\prime}_{1},\cos\theta^{\prime}_{2},\phi^{\prime}_{1},\phi^{\prime}_{2}$
as well as the angles in the center-of-mass (CM) frame of the final particles
$\cos\theta_{r}$ and $\phi_{r}$. Within each panel, we provide a comparison
between the distributions for phase-space (PS), NP, and PP channels.

In these distributions we fix the bins for one variable and integrate with
respect to all other variables (or sum over all remaining bins). As a
consequence the one-fold distributions on the different channels appears quite
similar. This is, in a way, an averaging effect on the hadronic tensor.
However, despite this overall similarity, the NP channel exhibits notable
distinctions from the NP and PP channels. The differences between the phase
space and the two emission channels are typically on the order of 10%–20%
depending on the kinematics.

We observe, when comparing the two top panels, that the distribution
$P(p^{\prime}_{1})$ is equal to the $P(p^{\prime}_{2})$ distribution for PP
emission. This equality arises due to the symmetry under the exchange of two
protons. The same symmetry holds for the phase space, but in this case, it’s
because, by definition, the phase-space distribution depends solely on
kinematics. However, in the case of NP emission, this symmetry is not observed
due to the inherent differences between the two outgoing particles and the
isospin dependence of the hadronic tensor. Recall that particle 1 is a neutron
and particle 2 is a proton. Consequently, the neutron distribution exhibits a
more pronounced peak at higher momenta compared to the proton distribution.

Moving on to the panels of the second row in Fig. 10, we examine the
distributions with respect to the polar emission angles,
$P(\cos\theta^{\prime}_{1})$ and $P(\cos\theta^{\prime}_{2})$. Ideally, these
two distributions should be exactly equal for the phase space, but in the
figure, they do not appear exactly identical one being slightly higher than
the other. This small discrepancy results from the numerical error introduced
when discretizing the integrals into bins and the calculation method not
treating particles 1 and 2 symmetrically. Consequently, the differences
between these curves provide an estimate of the numerical error incurred when
calculating these angular distributions, which here is approximately 5%. To
reduce this error, one would need to decrease the bin size accordingly.
However, this would significantly increase the number of exclusive events that
need to be summed (and computed). It’s important to keep this error in mind
when interpreting the results of this study. Looking at the probability
distributions in the panels of the second row, we notice that they increase
with $\cos\theta^{\prime}_{i}$, which means they decrease with the angle
itself. This indicates that forward emission (with forward being the direction
of the neutrino) is more probable, while backward emission is less likely.

Now we examine the two panels of the third row in Fig. 10 that display the
distributions over the azimuthal angles. This indicates the probability of the
orientation of the reaction plane of particle one with respect to the
scattering plane when the other emerges in any direction, or vice versa in the
case of particle 2. We observe that the prevailing trend is to be emitted
predominantly in the semi-plane with an angle of $\pi$, which is where the
momentum transfer vector is contained (see Fig. 2). There is an asymmetry in
the case of NP emission when switching particle 1 for particle 2. This results
in the neutron (particle 1) having a greater tendency to be emitted in the
$\pi$-plane (and adjacent planes) than the proton (particle 2).

Let’s now examine the distributions with respect to the emission angles in the
CM of the two final particles, shown in the bottom panels of Fig. 10. The
distribution of $\cos\theta_{r}$ is quite flat in the case of phase space and
PP emission, indicating that PP emission is approximately consistent with an
isotropic distribution in the center of mass. However, it’s also evident that
NP emission is not as compatible with this hypothesis, as the angular
distribution deviates from the phase space distribution. It is larger for
forward angles and smaller for other angles.

In the case of the distribution of the relative azimuthal angle $\phi_{r}$,
both the PP and NP channels deviate form a pure isotropic distribution in the
CM. The NP channel stands out as distinctly different from the other two. In
this case, both the distributions for the PS and PP channels exhibit two
maxima around $\pi/2$ and $3\pi/2$, with a minimum at $\pi$. In contrast, the
NP distribution displays two larger maxima that are closer to $\pi$.

Figure 11: The same of Fig. 10 for kinematics 2.

Similar differences are observed in Figure 11, which corresponds to kinematics
2. However, comparing with Figure 10, it becomes evident that the specific
shape of the distribution also whimsically depends on the leptonic kinematics.
Nonetheless, in almost all cases, the distributions in the case of NP emission
are notably different from the case of PP emission.

Figure 12:  One-fold distributions for semi-inclusive two-nucleon emission
induced by antineutrinos for kinematics 1\. In each panel, we compare the
distributions for phase-space (PS), NP, and NN emission.  Figure 13: The same
as Fig. 12 for kinematics 2.

To complete this discussion, we present in Figs. 12 and 13 results for
antineutrino scattering for kinematics 1 and 2, to be compared with the
corresponding figures 10 and 11 for neutrinos. They show similar trends, with
the most noticeable differences occurring again for NP emission as compared to
the NN emission and the phase space distribution, which are more similar. The
most notable differences between NP and the others are observed in the
distributions with respect to the azimuthal angles.

In any case, it is notable that the significant differences introduced by the
hadronic tensor in the two-fold distributions concerning relative angles, as
shown in Figures 4 and 5, tend to smooth out when we transition to
distributions of the observable variables related to the momenta of the two
nucleons. In other words, the very pronounced differences arising from
considering different emission channels compared to the isotropic distribution
in phase space seem to lessen when we consider the physical variables as a
result of an average over bins. It is expected that these differences become
much more pronounced when considering the complete six-fold distribution,
where the full impact of the hadronic tensor should be more prominent.

##  V Conclusions

In this work, we have undertaken a comprehensive study of the semi-inclusive
cross section for two-nucleon emission induced by neutrinos and antineutrinos
within a RMF model of nuclear matter, including MEC. The primary aim of this
study was to assess the differences between various emission channels, namely
NP, PP, and NN, and to compare them with the pure phase space model commonly
employed in neutrino Monte Carlo event generators.

The novelty in this work lies in the inclusion of the full semi-inclusive
hadronic tensor, derived from a microscopic calculation, whose integral
recovers the inclusive 2p2h hadronic tensor. Given the complexity of a
detailed study of the full six-dimensional distribution dependent on the final
nucleon momenta, involving numerous variables, we conducted this preliminary
investigation focusing on partial distributions dependent on only one or two
variables (one-fold and two-fold distributions). We integrated over the
remaining variables to gain insights into how the differences imposed by the
hadronic tensor propagate through these observables, for various lepton
kinematics.

Our findings reveal significant variations between the different emission
channels, especially when compared to the pure phase space model. These
differences are particularly pronounced in distributions involving the
relative angles in the CM frame. This is evidence that the isotropic
distribution in the center of mass undergoes significant modifications when
the hadronic tensor is included, effectively making it non-isotropic. This is
particularly relevant in the case of NP emission, which seems to be related to
the fact that the NP cross section is smaller than the PP cross section and is
thus more sensitive to the influence of the hadronic tensor.

In the case of two-fold distributions involving variables of the two emitted
particles, the differences with the pure phase-space model are smaller due to
the averaging effect of integrating over the remaining variables, which
smoothes out the differences. Furthermore, these distributions are largely
determined by kinematics, as energy and momentum conservation plays a
significant role, and the hadronic tensor mainly influences the finer details
of each distribution without altering the fundamental kinematic constraints.

From our results, we have observed strong correlations between pairs of
variables $(p^{\prime}_{1},p^{\prime}_{2})$,
$(\theta^{\prime}_{1},\theta^{\prime}_{2})$, and
$(\phi^{\prime}_{1},\phi^{\prime}_{2})$, which is evident from the
corresponding two-fold distributions. When considering other combinations of
two variables, there is generally very little correlation between them.

Appreciable differences are observed between the NP and PP distributions when
the hadronic tensor is included in neutrino scattering. One notable result is
that when an NP pair is emitted with neutrinos, it is more likely for the
neutron to carry more energy than the proton, whereas the opposite is true for
antineutrino scattering.

We have also presented results for the one-fold distributions, which were
obtained by fixing one semi-inclusive variable and integrating over the
remaining ones. These distributions undergo relatively few variations when the
hadronic tensor is included, primarily due to the smoothing effect of
integration. Nevertheless, significant differences are observed between the
cases of NP and PP emission, for neutrino scattering, and NP and NN emission,
for antineutrino scattering, depending on the kinematics.

In conclusion, this work has explored semi-inclusive two-nucleon emission
reactions induced by neutrinos, considering a microscopic treatment of the
corresponding hadronic tensor. We have emphasized the differences with respect
to the approach that assumes an isotropic distribution of the outgoing
nucleons in the center of mass. This research is expected to be valuable for
analyses based on Monte Carlo event generators and offers the potential for
improving the reconstruction of incident neutrino energies. Future work could
extend this study to the full six-dimensional distribution, providing even
deeper insights into the impact of the hadronic tensor on these processes.

## Acknowledgments

Work supported by: Grant PID2020-114767GB-I00 funded by
MCIN/AEI/10.13039/501100011033; FEDER/Junta de Andalucia-Consejeria de
Transformacion Economica, Industria, Conocimiento y
Universidades/A-FQM-390-UGR20; and Junta de Andalucia (Grant No. FQM-225).

## References

  * (1) H. Gallagher, G. Garvey and G. P. Zeller, Ann. Rev. Nucl. Part. Sci. 61, 355-378 (2011).
  * (2) J. G. Morfin, J. Nieves and J. T. Sobczyk, Adv. High Energy Phys. 2012, 934597 (2012).
  * (3) J. A. Formaggio and G. P. Zeller, Rev. Mod. Phys. 84, 1307-1341 (2012).
  * (4) L. Alvarez-Ruso, Y. Hayato, and J. Nieves, New J. Phys. 16, 075015 (2014).
  * (5) U. Mosel, Ann. Rev. Nuc. Part. Sci. 66, 171 (2016).
  * (6) M. Sajjad Athar, A. Fatima and S. K. Singh, Prog. Part. Nucl. Phys. 129, 104019 (2023).
  * (7) A. Aguilar-Arevalo et al. (MiniBooNE Collaboration), Phys. Rev. D 81, 092005 (2010).
  * (8) A. Aguilar-Arevalo et al. (MiniBooNE Collaboration), Phys. Rev. D 88, 032001 (2013).
  * (9) G.A. Fiorentini et al. (MINERvA Collaboration) Phys. Rev. Lett. 111, 022502 (2013).
  * (10) K. Abe et al., (T2K Collaboration), Phys. Rev. D 87, 092003 (2013).
  * (11) M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rec. C 80, 065501 (2009).
  * (12) J. Nieves, I. Ruiz Simo, and M.J. Vicente Vacas, Phys. Rev. C 83, 045501 (2011).
  * (13) J.E. Amaro, M.B. Barbaro, J.A. Caballero, T.W. Donnelly, and C.F. Williamson, Phys. Lett. B 696, 151 (2011).
  * (14) O. Lalakulich, K. Gallmeister, and U. Mosel, Phys. Rev. C 86, 014614 (2012); 90, 029902(E) (2014).
  * (15) T. Van Cuyck, N. Jachowicz, R. Gonzalez-Jimenez, M. Martini, V. Pandey, J. Ryckebusch, and N. Van Dessel, Phys. Rev. C 94, no.2, 024611 (2016).
  * (16) T. Van Cuyck, N. Jachowicz, R. Gonzalez-Jimenez, J. Ryckebusch, and N. Van Dessel, Phys. Rev. C 95, no.5, 054611 (2017).
  * (17) E. Berrueta Martinez, A. Mariano, and C. Barbero, Phys. Rev. C 103, no.1, 015503 (2021).
  * (18) A. Lovato, S. Gandolfi, J. Carlson, S. C. Pieper, and R. Schiavilla, Phys. Rev. C 91, no.6, 062501 (2015).
  * (19) A. Lovato, J. Carlson, S. Gandolfi, N. Rocco, and R. Schiavilla, Phys. Rev. X 10, no.3, 031068 (2020).
  * (20) N. Rocco, C. Barbieri, O. Benhar, A. De Pace, and A. Lovato, Phys. Rev. C 99, no.2, 025502 (2019).
  * (21) R. Gran, J. Nieves, F. Sanchez and M.J. Vicente Vacas, Phys. Rev. D 88, 113007 (2013).
  * (22) M. Martini and M. Ericson, Phys. Rev. C 87, 065501 (2013).
  * (23) M. Martini and M. Ericson, Phys. Rev. C 90, 025501 (2014).
  * (24) J.E. Amaro, M.B. Barbaro, J.A. Caballero and T.W. Donnelly, Phys. Rev. Lett. 108, 152501 (2012).
  * (25) J.M. Campbell, et. al., [arXiv:2203.11110 [hep-ph]].
  * (26) J.T. Sobczyk, Phys. Rev. C 86, 015504 (2012).
  * (27) C. Andreopoulos, A. Bell, D. Bhattacharya, F. Cavanna, J. Dobson, S. Dytman, H. Gallagher and P. Guzowski et al., Nucl. Instrum. Meth. A 614, 87 (2010).
  * (28) Teppei Katori, AIP Conf. Proc. 1663, 030001 (2015).
  * (29) S. Dolan, G. D. Megias and S. Bolognesi, Phys. Rev. D 101, no.3, 033003 (2020).
  * (30) S. Dytman, Y. Hayato, R. Raboanary, J. T. Sobczyk, J. Tena Vidal and N. Vololoniaina, Phys. Rev. D 104, no.5, 053006 (2021).
  * (31) J. M. Franco-Patino, R. Gonzalez-Jimenez, S. Dolan, M. B. Barbaro, J. A. Caballero, G. D. Megias and J. M. Udias, Phys. Rev. D 106, no.11, 113005 (2022).
  * (32) J. E. Sobczyk, J. Nieves and F. Sanchez, Phys. Rev. C 102, no.2, 024601 (2020).
  * (33) K. Abe et al. [T2K], Phys. Rev. D 98, no.3, 032003 (2018).
  * (34) X. G. Lu et al. [MINERvA], Phys. Rev. Lett. 121, no.2, 022504 (2018).
  * (35) T. Cai et al. [MINERvA], Phys. Rev. D 101, no.9, 092001 (2020).
  * (36) P. Abratenko et al. [MicroBooNE], Phys. Rev. D 102, no.11, 112013 (2020).
  * (37) P. Abratenko et al. [MicroBooNE], Phys. Rev. Lett. 125, no.20, 201803 (2020).
  * (38) O. Moreno, T. W. Donnelly, J. W. Van Orden and W. P. Ford, Phys. Rev. D 90, no.1, 013014 (2014).
  * (39) J. W. Van Orden and T. W. Donnelly, Phys. Rev. C 100, no.4, 044620 (2019).
  * (40) J. M. Franco-Patino, J. Gonzalez-Rosa, J. A. Caballero and M. B. Barbaro, Phys. Rev. C 102, no.6, 064626 (2020).
  * (41) J. M. Franco-Patino, M. B. Barbaro, J. A. Caballero and G. D. Megias, Phys. Rev. D 104, no.7, 073008 (2021).
  * (42) M. B. Barbaro, PoS NuFact2021, 007 (2022). [arXiv:2203.16132 [hep-ph]].
  * (43) J. M. Franco-Patino, M. B. Barbaro, J. A. Caballero and G. Megias, PoS NuFact2021, 227 (2022).
  * (44) J. M. Franco-Patino, S. Dolan, R. Gonzalez-Jimenez, M. B. Barbaro, J. A. Caballero and G. D. Megias, [arXiv:2304.01916 [hep-ex]].
  * (45) A. De Pace, M. Nardi, W. M. Alberico, T. W. Donnelly and A. Molinari, Nucl. Phys. A 726, 303 (2003).
  * (46) I. Ruiz Simo, J. E. Amaro, M. B. Barbaro, A. De Pace, J. A. Caballero and T. W. Donnelly, J. Phys. G 44, no.6, 065105 (2017).
  * (47) R. Rosenfelder, Ann. Phys. (N.Y.) 128, 188 (1980).
  * (48) B. D. Serot and J. D. Walecka, in _Advances in Nuclear Physics_ , edited by J. W. Negele and E. Vogt (Plenum, New York, 1986), Vol. 16.
  * (49) K. Wehrberger, Phys. Rep. 225, 273 (1993).
  * (50) G. D. Megias, J. E. Amaro, M. B. Barbaro, J. A. Caballero and T. W. Donnelly, Phys. Rev. D 94, no.1, 013012 (2016).
  * (51) G. D. Megias, J. E. Amaro, M. B. Barbaro, J. A. Caballero, T. W. Donnelly and I. Ruiz Simo, Phys. Rev. D 94, no.9, 093004 (2016).
  * (52) G. D. Megias, M. B. Barbaro, J. A. Caballero, J. E. Amaro, T. W. Donnelly, I. Ruiz Simo and J. W. Van Orden, J. Phys. G 46, no.1, 015104 (2019).
  * (53) M. V. Ivanov, A. N. Antonov, G. D. Megias, J. A. Caballero, M. B. Barbaro, J. E. Amaro, I. Ruiz Simo, T. W. Donnelly and J. M. Udias, Phys. Rev. C 99, no.1, 014610 (2019).
  * (54) V. L. Martinez-Consentino, I. R. Simo and J. E. Amaro, Phys. Rev. C 104, no.2, 025501 (2021).
  * (55) V. L. Martinez-Consentino, J. E. Amaro and I. Ruiz Simo, Phys. Rev. D 104, no.11, 113006 (2021).
  * (56) V. L. Martinez-Consentino and J. E. Amaro, [arXiv:2308.13847 [nucl-th]].
  * (57) P. R. Casale, J. E. Amaro and M. B. Barbaro, Symmetry 15 1709 (2023).
  * (58) I. Ruiz Simo, J. E. Amaro, M. B. Barbaro, A. De Pace, J. A. Caballero and T. W. Donnelly, J.Phys. G44, no.6, 065105 (2017).
  * (59) E. Hernandez, J. Nieves and M. Valverde, Phys. Rev. D 76, no.3, 033005 (2007)
  * (60) I. Ruiz Simo, C. Albertus, J. E. Amaro, M. B. Barbaro, J. A. Caballero and T. W. Donnelly, Phys. Rev. D 90, no.5, 053010 (2014).
  * (61) I. Ruiz Simo, J. E. Amaro, M. B. Barbaro, A. De Pace, J. A. Caballero, G. D. Megias and T. W. Donnelly, Phys. Lett. B 762, 124-130 (2016).
  * (62) V. L. Martinez-Consentino, J. E. Amaro, P. R. Casale and I. Ruiz Simo, Phys. Rev. D 108, no.1, 013007 (2023).