# ACHILLES: A novel event generator for electron- and neutrino-nucleus
scattering

Joshua Isaacson${}^{\,{\rm a}}$, William I. Jay${}^{\,{\rm d}}$, Alessandro
Lovato${}^{\,{\rm b,c}}$, Pedro A. N. Machado${}^{\,{\rm a}}$, Noemi
Rocco${}^{\,{\rm a}}$ ${}^{\,{\rm a}}$Theoretical Physics Department, Fermi
National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA
${}^{\,{\rm b}}$Physics Division, Argonne National Laboratory, Argonne,
Illinois 60439, USA
${}^{\,{\rm c}}$INFN-TIFPA Trento Institute of Fundamental Physics and
Applications, Via Sommarive, 14, 38123 Trento, Italy
${}^{\,{\rm d}}$Center for Theoretical Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA

###### Abstract

We present a novel lepton-nucleus event generator: ACHILLES, A CHIcagoLand
Lepton Event Simulator. The generator factorizes the primary interaction from
the propagation of hadrons in the nucleus, which allows for a great deal of
modularity, facilitating further improvements and interfaces with existing
codes. We validate our generator against high quality electron-carbon
scattering data in the quasielastic regime, including the recent CLAS/e4v
reanalysis of existing data. We find good agreement in both inclusive and
exclusive distributions. By varying the assumptions on the propagation of
knocked out nucleons throughout the nucleus, we estimate a component of
theoretical uncertainties. We also propose novel observables that will allow
for further testing of lepton-nucleus scattering models. ACHILLES is readily
extendable to generate neutrino-nucleus scattering events.

###### pacs:

24.10.Cn,25.30.Pt,26.60.-c

††preprint: FERMILAB-PUB-22-411-T, MIT-CTP/5428

##  I Introduction

Interactions between leptons with nuclei at beam energies in the 10 MeV$-$10
GeV range are key to properly interpret a multitude of experiments, e.g. those
that probe neutrino oscillations Itow  _et al._ (2001); Ayres  _et al._
(2004); Acciarri  _et al._ (2015); Abe  _et al._ (2018), electron-nucleus
scattering Abbott  _et al._ (1998); Dutta  _et al._ (2000); Mecking  _et al._
(2003); Dutta  _et al._ (2003); Rohe  _et al._ (2005); Burkert  _et al._
(2020), dark matter and dark sectors or even muon-specific new gauge forces
Battaglieri  _et al._ (2017); Kahn  _et al._ (2018); Åkesson  _et al._
(2018a); Battaglieri  _et al._ (2019); Agrawal  _et al._ (2021); Battaglieri
_et al._ (2021). Nevertheless, modeling these interactions with the percent-
level precision required by experimental analyses Acciarri  _et al._ (2015) is
a formidable challenge. These difficulties are primarily due to non-
perturbative nuclear dynamics, which play an important role in both the hard
interaction vertex and in the propagation of the struck nucleons before they
exit the nucleus Serber (1947); Metropolis  _et al._ (1958); Bertini (1963);
Cugnon (1980); Bertsch  _et al._ (1984); Stoecker and Greiner (1986); Bauer
_et al._ (1986); Bertsch and Das Gupta (1988); Botermans and Malfliet (1990);
Cassing  _et al._ (1990); Danielewicz and Bertsch (1991); Teis  _et al._
(1997); Cugnon  _et al._ (1997); Golubeva  _et al._ (1998); Boudard  _et al._
(2002); Hayato (2002); Casper (2002); Duarte (2007); Andreopoulos  _et al._
(2010); Iwamoto  _et al._ (2010); Buss  _et al._ (2012); Golan  _et al._
(2012); Sawada  _et al._ (2012); Uozumi  _et al._ (2012); Battistoni  _et al._
(2013, 2015); Mosel (2019); Niewczas and Sobczyk (2019); Isaacson  _et al._
(2021a); Dytman  _et al._ (2021).

Accurately modeling the hard-scattering cross section between a lepton and a
nuclear constituent presents non-trivial difficulties. The individual
couplings between leptons and nucleons are parametrized in terms of single-
nucleon form factors, which depend on the momentum transfer associated with
the process and involve non-perturbative QCD dynamics. These form factors can
be either calculated from first principles using lattice QCD or fitted to
experimental data (see e.g. Refs. Bhattacharya  _et al._ (2011); Meyer  _et
al._ (2016); Davoudi  _et al._ (2021); Borah  _et al._ (2020); Meyer  _et al._
(2022) and references therein). Whichever method is used, care must be taken
in properly estimating uncertainties, especially in the axial sector where
experimental data are scarce. In addition, scattering does not occur on a
collection of free nucleons (which could be described trivially within a
Fermi-gas model). Instead, the real-world target nucleus is a correlated
quantum many-body system; hence, a percent-level theoretical description of
scattering off of a bound nucleon must capture many-body correlation effects
within a well-defined factorization scheme like the impulse approximation and
its generalizations (see e.g. Refs. Benhar  _et al._ (2008); Martini  _et al._
(2009, 2010); Amaro  _et al._ (2011); Gran  _et al._ (2013); Benhar  _et al._
(2015)).

The final-state interactions (FSI) that the struck nucleons undergo before
exiting the nucleus are also extremely complex phenomena Serber (1947);
Metropolis  _et al._ (1958); Bertini (1963); Cugnon (1980), subject to non-
perturbative single and many-nucleon effects. Examples of the former are
nucleon excitations into a $\Delta$ isobar, which quickly decays into a pion-
nucleon state. The latter include correlations induced by realistic two- and
three-nucleon forces. A fully quantum mechanical description of these
processes presents an exponentially hard computational problem. Sophisticated
nuclear many-body methods leverage leadership-class computing resources to
tackle this real-time nuclear dynamics problem but are limited to inclusive
processes and to the non-relativistic regime Lovato  _et al._ (2020); Sobczyk
_et al._ (2021). Exploratory calculations on quantum devices show promise
Roggero and Carlson (2019); Roggero  _et al._ (2020) for treating fully-
exclusive processes. However, the inclusion of relativistic effects poses non-
trivial challenges, and their application to realistic systems seems to remain
a distant goal.

Over the years, a number of complementary methods have been developed to
capture the leading effects of FSI. The most sophisticated ones start from the
Kadanoff-Baym integro-differential equations for the evolution of the entire
nuclear system. In practice, state-of-the-art transport codes solve truncated
versions of these equations, which in turn make a proper estimation of theory
uncertainties much harder Kadanoff and Baym (1962); Botermans and Malfliet
(1990); Cassing  _et al._ (1990); Teis  _et al._ (1997); Buss  _et al._
(2012); Mosel (2019). On the other hand, intranuclear cascade (INC) approaches
approximately solve the transport equation by evaluating the collision term
stochastically. The main approximation of INC models is that of classical
propagation between consecutive quantum-mechanical scatterings. Hence, they
are applicable in the regime in which the de Broglie wavelength of the
nucleons is much smaller than the range of the interaction, which is in turn
smaller than the average distance between nucleons Cugnon and Henrotte (2003).
Hence, the applicability of INCs is in principle restricted to nucleons with
kinetic energies above $\sim 200$ MeV, although many observables in heavy-ion
reactions at less than $100$ MeV are reproduced well in practice Dore  _et
al._ (2001); Uozumi  _et al._ (2012).

Besides the aforementioned intrinsic limitations, INCs involve a number of
additional model-specific prescriptions. For instance, the kinematic variables
of a struck nucleon are drawn from a model of the nuclear ground-state;
typically the (local or global) Fermi gas or realistic spectral functions. In
addition, some assumptions for the propagation in the nuclear medium are made.
Examples of the latter are the use of mean free path estimates obtained from
the nuclear density or propagating nucleons as if they were “hard-spheres”
with a radius that is proportional to the square root of the total nucleon-
nucleon cross section. Finally, corrections due to the nuclear environment are
also typically included by means of average nuclear potentials and imposing
Pauli blocking in nucleon-nucleon collisions.

To assess the reliability of such _effective_ descriptions, INCs should be
systematically and extensively validated against available experimental data.
Electron-nuclei scattering experiments offer large, high-quality data samples
including a broad range of experimental observables, which can be used both to
benchmark different INC models and to gauge the reliability of their
assumptions. Due to the interplay of all these effects, modeling lepton-
nucleus interactions constitutes a remarkable challenge. Nevertheless, state-
of-the-art neutrino event generators have only recently started comparing
their predictions to electron scattering data Buss  _et al._ (2008); Leitner
_et al._ (2009); Mosel and Gallmeister (2019); Isaacson  _et al._ (2021a);
Ankowski and Friedland (2020); Dytman  _et al._ (2021); Khachatryan  _et al._
(2021). The results of these initial comparisons reveal that existing event
generators do not describe electron-nucleus scattering data, and consequently
neutrino-nucleus scattering data, to the precision level required by next-
generation experiments, such as DUNE Abi  _et al._ (2020a, b).

In this paper we take a first step towards a full-fledged lepton-nucleus event
generator: ACHILLES, A CHIcagoLand Lepton Event Simulator. Three aspects of
ACHILLES’ design are worth highlighting explicitly. First, as neutrino physics
enters the precision era, we expect that event generators will need to
incorporate many technical improvements and new physical insights. The need
for robust, quickly extensible codebases is therefore acute, presenting a
challenge not only in scientific computing but also in software engineering.
To help solve this inherent difficulty, one of the core design principles of
our code is _modularity_. We have endeavored to divide the code clearly into
individual pieces that describe the different physics processes within lepton-
nucleus scattering. Examples of the constituent parts include the description
of the initial state of the nucleus; the “hard scattering” between the lepton
and constituents of the nucleus; the intranuclear cascade process; and the
nuclear potential. The advantage of modularity is obvious: an improvement on a
specific part of the code requires minimal effort, since all parts are
independent. Well-implemented modular design thus facilitates future
improvements and keeping abreast with new advances in the description of
lepton-nucleus scattering. We have also tried to provide a user-friendly
interface to available tools for physics beyond the Standard Model (BSM), such
as the recent lepton-tensor interface developed by some of the current authors
Isaacson  _et al._ (2021b). The ultimate goal is to enable ACHILLES users to
proceed seamlessly from writing down a BSM Lagrangian to generating events.

Second, the physical structure of the scattering problem also imposes
important constraints. Neutrino-nucleus scattering involves both vector and
axial form factors (associated with one- and two-body nuclear electroweak
current operators), while electron-nucleus scattering is dominated by vector
form factors from photon exchange. Thus, any neutrino event generator should
be able to describe electron-nucleus scattering data as a special case of the
more general problem. Therefore, the benchmarking foundation for any lepton-
nucleus event generator must be extensive validation against inclusive and
semi-inclusive electron-nucleus scattering data. The present paper is an
attempt to begin laying this foundation for ACHILLES.

Third, to leverage the existing analysis tools developed and maintained by the
high energy event generators at the LHC, we adopt the HepMC3 output format
Buckley  _et al._ (2021). This format allows easy interface with analysis
tools such as Rivet Buckley  _et al._ (2013); Bierlich  _et al._ (2020) and
eventually Nuisance Stowell  _et al._ (2017), saving valuable research time
for users to focus more on physics and less on coding. The HepMC3 output
format permits arbitrary parameters to be added to an event, allowing
additional event information required by neutrino experiments to be included
in a simple and straightforward manner.

This paper constitutes the first step towards the full development of
ACHILLES. We perform a comparison between our model of electron-nuclei
interactions, including the INC model developed in Ref. Isaacson  _et al._
(2021a), against the recent reanalysis of electron-carbon scattering data by
the CLAS/e4v collaboration Khachatryan  _et al._ (2021). Compared to our
previous work, and to obtain more realistic results for exclusive observables,
we implement a nuclear potential and simulate the propagation of nucleons
within this potential. We focus on the inclusive quasielastic (QE) cross
section, which is better understood than other cross section channels at the
energies of interest Rocco  _et al._ (2019); Ruso  _et al._ (2022). We also
consider the angular dependent proton yield, as well as a few other
kinematical observables in the QE regime.

The comparisons performed here test the modelling of the cross section; the
impact of initial-state nuclear configurations, particularly those obtained
via quantum Monte Carlo methods Carlson  _et al._ (2015) and the roles of the
spectral function, in-medium modifications, and Pauli blocking effects. Taken
together, these comparisons provide valuable insight to the physics of
intranuclear cascades. Besides these comparisons, we also propose new
observables that may help in further testing models of electron-nuclei and
neutrino-nuclei interactions. All proposed observables can be readily
extracted from existing CLAS data.

The paper is organized as follows. Sec. II lays out general considerations for
lepton-nucleus scattering. Sec. III compares the results of ACHILLES to
experimental data on inclusive electron-nucleus scattering. Sec. IV.1
describes in medium effects from the nuclear potential. Sec. V provides
comparisons to exclusive data. Sec. VI proposes novel observables to further
test the interaction modelling, following by conclusions in Sec. VII.

##  II General lepton-nucleus scattering

The general expression of the differential cross section for a scattering
process involving a target nucleus $A$ and a lepton $\ell$ leading to a given
final state reads

$\displaystyle d\sigma=$
$\displaystyle\left(\frac{1}{|v_{A}-v_{\ell}|}\frac{1}{4E_{A}^{\rm
in}E_{\ell}^{\rm in}}\right)\big{|}\mathcal{M}\big{|}^{2}$
$\displaystyle\times\prod_{f}\frac{d^{3}p_{f}}{(2\pi)^{3}}(2\pi)^{4}\delta^{4}\Big{(}k_{A}+k_{\ell}-\sum_{f}p_{f}\Big{)}\,.$
(1)

We denote the initial- and final-state momenta by

$\displaystyle k_{\ell}^{\mu}=(E_{\ell}^{\rm in},{\bf
k})\quad\text{incoming}~{}\ell,$ (2) $\displaystyle k_{A}^{\mu}=(E_{A}^{\rm
in},{\bf k}_{A})\quad\text{incoming}~{}A,$ (3) $\displaystyle
p_{f}^{\mu}=(E_{f},{\bf p}_{f})\quad\text{outgoing~{}particle~{}}f,$ (4)

where the index $f$ refers to all the possible hadronic and leptonic final
state particles. The first term in parenthesis is the flux of incoming
particles, with $v_{\ell,A}$ being the velocities, the second encodes the
matrix element, and the last line is the phase space for the outgoing
particles. In the one-boson exchange approximation, the squared amplitude
reads

$\big{|}\mathcal{M}\big{|}^{2}=L_{\mu\nu}W^{\mu\nu}\frac{1}{P^{2}},$ (5)

where $P$ is a generic vector boson propagator, while $L$ and $W$ denote the
leptonic and hadronic tensors, respectively. The leptonic tensor is completely
determined by the leptonic process (e.g. neutrino charged-current
interactions). The hadronic tensor on the other hand contains all information
on nuclear dynamics and it is expressed as

$W_{\mu\nu}=\langle\Psi_{0}|J^{\dagger}_{\mu}(q)|\Psi_{f}\rangle\langle\Psi_{f}|J_{\nu}(q)|\Psi_{0}\rangle\,,$
(6)

where $\Psi_{0}$ and $\Psi_{f}$ denotes the hadronic initial and final states,
respectively.

Carrying out the full calculation of the many-body wave function and its real-
time evolution is an exponentially-hard computational problem. To handle it,
the reaction process is modeled by separating the primary interaction vertex
from the propagation of the struck particles out of the nucleus.
Schematically, this division can be expressed by considering the full matrix
element squared as:

$\displaystyle|\mathcal{M}(\\{k\\}\rightarrow\\{p\\})|^{2}=\big{|}\sum\int_{\;\;p^{\prime}}$
$\displaystyle\mathcal{V}(\\{k\\}\rightarrow\\{p^{\prime}\\})$
$\displaystyle\times\mathcal{P}(\\{p^{\prime}\\}\rightarrow\\{p\\})\big{|}^{2},$
(7)

where the $\\{k\\}(\\{p\\})$ is the set of all initial(final) state particle
momenta, $\mathcal{V}$ represents the primary interaction vertex producing
intermediate particles with momenta $\\{p^{\prime}\\}$, and $\mathcal{P}$
denotes the time evolution of the intermediate states to the final states
outside the nucleus. Calculating this equation exactly requires retaining full
quantum mechanical interference between the primary interaction vertex and the
subsequent re-interactions. Traditionally, due to the complexity of solving
Eq. (II) exactly, the calculation factorizes the two-step process as an
incoherent product:

$\displaystyle|\mathcal{M}(\\{k\\}\rightarrow\\{p\\})|^{2}\simeq\sum\int_{\;\;p^{\prime}}$
$\displaystyle|\mathcal{V}(\\{k\\}\rightarrow\\{p^{\prime}\\})|^{2}$
$\displaystyle\times|\mathcal{P}(\\{p^{\prime}\\}\rightarrow\\{p\\})|^{2}.$
(8)

This treatment is similar to the approach taken by the collider community when
dressing hard-scattering cross sections with parton showers (see e.g. Ref.
Höche (2015) for a pedagogical discussion). By construction, this
approximation neglects inference between primary interaction vertices that
give rise to identical final states while leaving inclusive observables
unaffected. It is expected that these interference effects are subdominant,
and a detailed investigation is left to a future work. In ACHILLES, the
subsequent evolution probability $|\mathcal{P}|^{2}$ is handled semi-
classically using the algorithm developed in Ref. Isaacson  _et al._ (2021a).

Eqs. (II) and (II) retain full generality for lepton-nucleus scattering, but
implementing them in a concrete calculation requires several choices about the
relevant degrees of freedom. First, one must specify the initial-state nuclear
constituents $\\{k\\}$ which participate in the vertex $\mathcal{V}$. This
question is closely related to the choice of a factorization scheme, to which
the following section is dedicated. Second, one must specify the intermediate-
state particles which can appear, either from production at the primary
interaction vertex $\mathcal{V}$ or in the system’s subsequent evolution
$\mathcal{P}$. Briefly stated, the present work restricts to processes in
which protons and neutrons are the only active degrees of freedom. This choice
explicitly neglects, e.g., pion production at the primary interaction vertex.
Under this ansatz, the electroweak current of Eq. (6) is expanded as a sum of
one and two-nucleon operators

$J^{\mu}(q)=\sum_{i}j^{\mu}_{i}(q)+\sum_{i<j}j^{\mu}_{ij}(q)\,.$ (9)

Three- and higher-body terms have been found to be small Marcucci  _et al._
(2005) and are thus neglected here. Generalizing these expressions to other
mediators, such as scalars, is straightforward.

###  II.1 Factorization scheme

As mentioned above, this work focuses on a lepton scattering on a nucleus in
the quasielastic regime, in which the dominant reaction mechanism is assumed
to be single-nucleon knockout: $\ell+A\rightarrow\ell^{\prime}+X$ where $\ell$
and $\ell^{\prime}$ denote the initial and final lepton states, $A$ is the
target nucleus, $X$ is the hadronic final state, which for example can be
composed of a single emitted nucleon and the remnant nucleus. We can rewrite
the cross section of Eq. (1) for a $2\rightarrow 2$ scattering process in a
background as

$\displaystyle d\sigma=$
$\displaystyle\frac{d^{3}p_{\ell}}{(2\pi)^{3}}\frac{1}{2E_{\ell}}\frac{d^{3}p_{X}}{(2\pi)^{3}}\frac{1}{2E_{X}}L_{\mu\nu}W^{\mu\nu}\frac{1}{P^{2}}$
(10)

where again $p$ refer to final state momenta, with $X$ denoting the hadronic
final state, and the factor $1/|v_{A}-v_{\ell}|=1$ in the lab frame
(neglecting the lepton mass). Note that, for convenience, we absorbed the
factors $1/2E_{A(\ell)}^{\rm in}$ in the hadronic (leptonic) tensor, and we
are also embedding the delta function in the hadronic tensor. For large enough
values of the momentum transfer, the virtual boson primarily interacts with
individual bound nucleons, so that the hadronic final state can be
approximated by the factorized expression

$|\Psi_{f}\rangle=|p\rangle\otimes|\Psi_{f}^{A-1}\rangle\,,$ (11)

where $|p\rangle$ is a plane wave describing the propagation of the final
state nucleon with momentum $|{\bf p}|$, while $|\Psi_{f}\rangle_{A-1}$
denotes the $(A-1)$-body spectator system, which can be either in a bound or
unbound state. In addition, since we are focusing on the primary interaction
vertex, we have dropped the prime in the intermediate variables.

Retaining the one-body current contribution only in Eq. (9), the incoherent
contribution to the hadronic tensor is given by

$\displaystyle W^{\mu\nu}({\bf q},\omega)$
$\displaystyle=\\!\\!\sum_{h\in\\{p,n\\}}\int\frac{d^{3}k_{h}}{(2\pi)^{3}}\frac{1}{2E_{h}^{\rm
in}}dE^{\prime}S_{h}({\bf k}_{h},E^{\prime})$ $\displaystyle\times\langle
k|{j_{1b}^{\mu}}^{\dagger}|p\rangle\langle
p|j_{1b}^{\nu}|k\rangle(2\pi)^{4}\delta^{3}({\bf k}_{h}+{\bf q}-{\bf p}_{h})$
$\displaystyle\times\delta(\omega-E^{\prime}+m_{N}-E_{h})\,$ (12)

In the previous equation, ${\bf q}={\bf p}_{h}-{\bf k}_{h}$ is the momentum
transfer, $k_{h}=(E_{h}^{\rm in},{\bf k}_{h})$ and $p_{h}=(E_{h},{\bf p}_{h})$
denote the energy and momentum of the initial and final nucleons,
respectively. Note that, for brevity, we have suppressed the subscript $h$ in
the bras and kets. The spectral function yields the probability distribution
of removing a “hole” nucleon $h\in\\{p,n\\}$ with momentum ${\bf k}_{h}$ from
the target nucleus, leaving the residual $(A-1)$ system with an excitation
energy $E^{\prime}$, and it is defined as Benhar and Pandharipande (1993)

$\displaystyle S_{h}($ $\displaystyle{\bf k}_{h},E^{\prime})=$ (13)
$\displaystyle\sum_{f_{A-1}}|\langle\Psi_{0}||k\rangle\otimes|\Psi_{f}^{A-1}\rangle|^{2}\delta(E^{\prime}+E_{0}^{A}-E_{f}^{A-1}),$

where the sum runs through the possible final states of the $A-1$ spectator
nucleons, which can either be bound or in the continuum.

The spectral function of finite nuclei is generally expressed as a sum of a
mean-field and a correlation contribution. The first one describes the low
momentum and removal-energy region, and is associated with the residual $A-1$
system being in a bound state. The correlation contribution includes unbound
$A-1$ states for the spectator system, in which at least one of the spectator
nucleons is in the continuum, and it provides strength in the high momentum
and energy region. The nuclear spectral function has been evaluated within
different semi-phenomenological Benhar  _et al._ (1994); Ivanov  _et al._
(2019) and ab-initio many-body methods Rocco and Barbieri (2018); Barbieri
_et al._ (2019), including quantum Monte Carlo Andreoli  _et al._ (2022).

The one employed in this work has been obtained within the correlated basis
function theory of Ref. Benhar  _et al._ (1994). The low momentum and energy
contribution is determined by adjusting mean-field calculations to reproduce
($e$, $e^{\prime}$p) scattering measurements. The correlation part is derived
within the Local Density Approximation by convoluting the correlation
component of the spectral function obtained within the correlated basis
function theory for isospin-symmetric nuclear matter for a given value of the
density. Additional details regarding the spectral function and, in
particular, corrections to the impulse approximation stemming from final-state
interactions appear in Sec. III.2.

The spectral function is normalized as

$\displaystyle\int\frac{d^{3}k_{h}}{(2\pi)^{3}}dE^{\prime}S_{h}({\bf
k}_{h},E^{\prime})=\begin{cases}Z,\quad\quad\,h={\rm p},\\\ A-Z,\,h={\rm
n}.\end{cases}$ (14)

where $Z$ denotes the number of protons in the nucleus. After applying the
factorization ansatz to the hadronic final state, the phase space factor of
Eq. (10) can be rewritten as

$\frac{d^{3}p_{f}}{(2\pi)^{3}}\frac{1}{2E_{f}}\rightarrow\frac{d^{3}p_{h}}{(2\pi)^{3}}\frac{1}{2E_{h}}\sum_{f_{A-1}}$
(15)

where the discrete sum over the states of the remnant nucleus is embedded in
the spectral function as shown in Eq. (13).

A complete estimate of the theoretical uncertainty associated with the cross
section calculation would require assessing the error in the many-body
calculation of the spectral function, the inputs used to describe the
interaction vertex (i.e. couplings, form factors), and the factorization of
the hadronic final state. Achieving this goal is highly nontrivial and has not
been included in this work but future developments are discussed in Sec. VII.

##  III Inclusive electron-nucleus scattering

###  III.1 Theoretical preliminaries

Now we proceed to the concrete calculation of the electron-nucleus cross
section, which will be the basis of all comparisons between ACHILLES and
electron-carbon scattering data. We focus first on comparisons between our
theoretical predictions and experimental data for the quasielastic inclusive
electron-12C cross section using the aforementioned factorization scheme and
spectral function formalism. For those kinematics in which FSI are expected to
be negligible, the inclusive cross section provides a benchmarking test for
the model of the primary interaction ($|\mathcal{V}|^{2}$ in Eq. II), since
this observable is unaffected by the semi-classical propagation in the nuclear
medium, and hence by the INC.

The inclusive double differential cross section for the scattering of an
electron on an at-rest nucleus via one-photon exchange is written as (see Eqs.
(2-4) for notation)

$\displaystyle\Big{(}\frac{d^{2}\sigma}{dE_{e}d\Omega_{e}}\Big{)}$
$\displaystyle=\frac{\alpha^{2}}{Q^{4}}{E_{e}^{2}}L_{\mu\nu}W^{\mu\nu}\,,$
(16)

where $\alpha\simeq 1/137$ is the fine structure constant, and $\Omega_{e}$ is
the scattering solid angle in the direction specified by ${\bf p}_{e}$. The
energy and the momentum transfer are denoted by $\omega$ and q, respectively,
with $Q^{2}=-q^{2}={\bf q}^{2}-\omega^{2}$. The lepton tensor is fully
determined by the lepton kinematic variables and, neglecting the electron
mass, it is given by

$L^{\mu\nu}=\frac{1}{E_{e}^{\rm
in}E_{e}}(k_{e}^{\mu}p_{e}^{\nu}+p_{e}^{\mu}k_{e}^{\nu}-g^{\mu\nu}\,k_{e}\cdot
p_{e})\,.$ (17)

The one-body electromagnetic current operator entering Eq. (12) is written as

$\displaystyle
j^{\mu}_{1b}={\mathcal{F}}_{1}\gamma^{\mu}+i\sigma^{\mu\nu}q_{\nu}\frac{{\mathcal{F}}_{2}}{2m_{N}}\
,$ (18)

where the isoscalar (S) and isovector (V) form factors, ${\mathcal{F}}_{1}$
and ${\mathcal{F}}_{2}$, are given by combination of the Dirac and Pauli ones,
$F_{1}$ and $F_{2}$, as

$\displaystyle{\mathcal{F}}_{1,2}=\frac{1}{2}[F_{1,2}^{S}+F_{1,2}^{V}\tau_{z}]\,,$
(19)

$\tau_{z}$ is the isospin operator, and

$\displaystyle F_{1,2}^{S}=F_{1,2}^{p}+F_{1,2}^{n},\ \ \ \ \ \ \ \
F_{1,2}^{V}=F_{1,2}^{p}-F_{1,2}^{n}\,.$ (20)

The Dirac and Pauli form factors can be expressed in terms of the electric and
magnetic form factors of the proton and neutron as

$\displaystyle F_{1}^{p,n}=$ $\displaystyle\frac{G_{E}^{p,n}+\tau
G_{M}^{p,n}}{1+\tau},\ \ \ \ \ \ \ \
F_{2}^{p,n}=\frac{G_{M}^{p,n}-G_{E}^{p,n}}{1+\tau}$ (21)

with $\tau=Q^{2}/4m_{N}^{2}$. Therefore, the electromagnetic current can be
schematically written as $j^{\mu}_{1b,\rm
EM}=j^{\mu}_{\gamma,S}+j^{\mu}_{\gamma,z}$ where the first is the isoscar term
and the second is the isovector multiplied by the isospin operators
$\tau_{z}$. The above set of equations can be readily extended to the
electroweak case and higher multiplicity processes; an automation for
arbitrary leptonic tensors was developed in Isaacson  _et al._ (2021b).

The use of a realistic spectral function combined with a factorization scheme
has proven to reproduce a large fraction of the available electron scattering
data (see Ref. Rocco (2020) and references therein). Over the past few years,
the factorization scheme has been extended to account for two-nucleon currents
and pion-production mechanisms Benhar  _et al._ (2015); Ruiz Simo  _et al._
(2017); Hernandez  _et al._ (2007); Kamano  _et al._ (2013); Nakamura  _et
al._ (2015); Kamano  _et al._ (2016). The focus of the present work is the
quasielastic region, and we leave the implementation of additional channels to
a future work.

###  III.2 Comparison to data

The first comparison between ACHILLES and data can be found in Fig. 1 (for
technical details of ACHILLES, see App. A). We present the ACHILLES inclusive
$e$-C quasielastic cross section (red histogram) against data as a function of
the energy transfer $\omega$. Data is taken from several experiments at
different incoming electron energy and outgoing electron angle, from top left
to bottom right: 730 MeV and 37∘ O’Connell  _et al._ (1987); 961 MeV and 37.5∘
Sealock  _et al._ (1989); 1300 MeVand 37.5∘ Sealock  _et al._ (1989); 2500 MeV
and 15∘ Zeller (1973). In all four cases, the first peak, which is dominated
by quasielastic scattering, is quite well described by ACHILLES. Note that
meson-exchange currents provide additional strength in the dip region between
the quasielastic and the resonance peak Rocco  _et al._ (2019). The second
peak has large contributions from resonance production, a mechanism which has
not yet been implemented in ACHILLES, and therefore it is not expected to be
reproduced by the present version of the code.

Figure 1: Comparison of the ACHILLES event generator to electron-carbon
scattering data. Top left: Scattering with an incoming energy of 730 MeV at an
angle of 37∘, data is from O’Connell  _et al._ (1987). Top right: Scattering
with an incoming energy of 961 MeV at an angle of 37.5∘, data is from Sealock
_et al._ (1989). Bottom left: Scattering with an incoming energy of 1300 MeV
at an angle of 37.5∘, data is from Sealock  _et al._ (1989). Bottom right:
Scattering with an incoming energy of 2500 MeV at an angle of 15∘, data is
from Zeller (1973).

Given the large $Q^{2}$ values of the data displayed in Fig. 1, FSI between
the struck nucleon and the remnant nucleus are expected to be small and have
been neglected in the initial hard interaction. For kinematics in which the
factorization scheme is not expected to hold, different approaches have been
developed to account for quantum-mechanical effects in FSI in the quasi
elastic region. To correct the factorization scheme and spectral function
results, the real part of an optical nuclear potential Cooper  _et al._ (2009)
is added to the free energy spectrum of the outgoing nucleon and the cross
section is convoluted with a _folding function_ to account for rescattering
effects Benhar (2013); Benhar  _et al._ (2008). In the Relativistic Mean Field
approach, FSI between the outgoing nucleon and the residual nucleus are
accounted for by solving the associated Dirac equation using the same mean
field as used for the bound nucleon González-Jiménez  _et al._ (2019).
Including these corrections modifies the inclusive cross section, shifting the
quasi elastic peak to lower energy transfers and redistributing the cross
section strength to the high-energy-transfer tail Ankowski  _et al._ (2015).

However, the aforementioned approaches do not allow for an accurate treatment
of exclusive processes. In this regard, INCs are a common tool Boudard  _et
al._ (2002); Hayato (2002); Casper (2002); Andreopoulos  _et al._ (2010);
Golan  _et al._ (2012); Battistoni  _et al._ (2013, 2015); Isaacson  _et al._
(2021a); Dytman  _et al._ (2021); Ershova  _et al._ (2022) for modeling the
total hadronic state that escapes the nucleus after the hard interaction
vertex, as described by Eq. (II), and that could be observed in the detector.
INCs use probabilities to determine if an additional scattering occurs. This
probabilistic treatment, at least in current algorithms, neglects interference
effects. Therefore, by definition INCs leave inclusive observables, such as
the differential cross section displayed in Fig. 1, unchanged.

It is important to note that the FSI modeled by folding functions and the FSI
modeled by INCs arise from the same physics. The major differences between the
two approaches are the approximations used to include the imaginary part of
the nuclear potential. Folding functions account for the effects of FSI
including interference effects at the cost of integrating out information on
the final state nucleons. On the other hand, INCs capture the exclusive final
state nucleons at the cost of neglecting the interference effects. Therefore,
combining the two calculations in a single code results in effectively double
counting the imaginary part of the nuclear potential. Implementing
interference effects into INCs is beyond the scope of this work.

Having established that our interaction model of quasielastic interactions
(i.e., $\mathcal{V}$ in Eq. II) reproduces the experimental data, we next move
to comparisons with exclusive observables.

##  IV Intranuclear Cascade

Simulating the propagation of the nucleon involved in the hard scattering out
of the nucleus is a vital component of a neutrino event generator.
Intranuclear cascade models are a class of algorithms used to reproduce the
imaginary part of the nuclear potential using stochastic Monte Carlo methods.
Traditional techniques do not capture the quantum mechanical components
involved in this process. Recently, a new technique for intranuclear cascades
has been proposed in Ref. Isaacson  _et al._ (2021a) to begin capturing these
effects. Figure 2 shows the programmatic flow of our cascade model.

In this algorithm, the spatial distribution of neutron and protons are sampled
from nuclear configurations obtained from QMC calculations fully retaining
correlations effects. Their initial momentum is generated according to a local
Fermi gas model. Once the target and the projectile are initialized, the
particles are propagated using relativistic kinematics. In the simplest
approximation, these particles follow straight-lines trajectories, but an
option to bend these trajectories using nuclear potentials is discussed in the
next section. We follow a time-like approach for the propagation; at each step
of the propagation $\delta t$ we check if an interaction occurred according to
the nucleon-nucleon scattering cross section using either a Gaussian or
cylindrical probability model depending on the impact parameter. Originally,
the only in-medium effect was taken to be the Pauli principle, and below we
discuss updates using the nuclear potential. We keep two separate lists of
“propagating” and “spectators” particles. At the beginning of the event, the
projectile is the only propagating particle. Afterwards, each particle that
has collided with a spectator is promoted to a propagating one, while all the
others are still labeled as spectators. The particles are propagated until
they reach the surface of the nucleus where the nucleon is either recaptured
or escapes, based on its energy.

Figure 2:  The proposed algorithm for the INC model. This figure is reproduced
from Ref. Isaacson  _et al._ (2021a).

###  IV.1 Nuclear potential

In electron-nucleus and neutrino-nucleus scatterings, inclusive quantities may
be well described without detailed modelling of what happens when nucleons are
propagating out of the nucleus. The description of exclusive quantities is
more demanding. While nucleon-nucleon interactions are possibly the most
important effects to include in an INC model, the presence of a mean-field
nuclear potential may trap struck nucleons or deflect their trajectory,
effectively changing the number, momentum and direction of outgoing particles.
To account for this effect, we have implemented two different options as a
background potential, which depends on both the position $r$ and momentum $p$
of the propagating nucleon. Note that, in our approaches, only the real part
of the potential is included, since the imaginary part is captured by the hard
scattering in the intranuclear cascade. A similar approach of including the
potential into cascades was studied in Ref. Nikolakopoulos  _et al._ (2022).
The first potential considered is a non-relativistic potential defined by a
three-parameter fit to single-particle energy of infinite nuclear matter
Wiringa (1988), which is consistent with the variational ground-state
calculations of Wiringa, Fiks, and Fabrocini (WFF) Wiringa  _et al._ (1988).
Its functional form is given as

$U(p^{\prime},r)=\alpha[\rho(r)]+\frac{\beta[\rho(r)]}{1+\left(p^{\prime}/\Lambda[\rho(r)]\right)^{2}},$
(22)

where $p^{\prime}$ is the modulus of the three momentum of the propagating
nucleon, while $\alpha,\beta$, and $\Lambda$ are fit to reproduce the single-
particle energy of nuclear matter as obtained from the Urbana $v_{14}$ \+ TNI
Hamiltonian, and $\rho(r)$ is the local nuclear density at radius $r$. The
values of the aforementioned variables are

$\displaystyle\alpha(\rho)$
$\displaystyle=15.52(\rho/\rho_{0})+24.93(\rho/\rho_{0})^{2}~{}\text{MeV},$
(23) $\displaystyle\beta(\rho)$
$\displaystyle=-116(\rho/\rho_{0})~{}\text{MeV},$ (24)
$\displaystyle\Lambda(\rho)$
$\displaystyle=3.29-0.373(\rho/\rho_{0})~{}\text{fm}^{-1},$ (25)

where $\rho_{0}=0.16~{}\text{fm}^{-3}$ is the saturation density of nuclear
matter.

The other potential we adopted is based on the work of Ref. Cooper  _et al._
(2009) where proton-nucleus elastic and reaction cross section data are fitted
to determine global proton-nucleus optical potentials for energies between 20
and 1040 MeV for several nuclear targets, including carbon. The fitting can be
done with potentials in a Dirac equation or Schroedinger equation. For the
former case, the Dirac equation was used in the form

$\\{{\bm{\alpha}}\cdot{\bf
p}^{\prime}+\beta[m+U_{s}(r,E^{\prime})]+U_{0}(r,E^{\prime})+V_{c}(r)]\\}\psi({\bf
r})=E^{\prime}\psi({\bf r})$

where $V_{c}(r)$ denotes the Coulomb potential at a given nuclear radius $r$,
which is either computed from Woods-Saxon-like charge distribution Clark  _et
al._ (2006) or taken from data when they are available, and $E^{\prime}$ is
the energy of the propagating nucleon. The quantities determined by the
fitting procedure are $U_{s}(r,E^{\prime})$ and $U_{0}(r,E^{\prime})$, the
scalar and vector optical potentials, respectively; they include a real and an
imaginary part. To obtain an effective optical potential for the Schroedinger
equation, it is helpful to write down a standard reduction of the Dirac
equation to second order form. The equation for the upper two components is

$\displaystyle\\{p^{\prime 2}+$ $\displaystyle 2E^{\prime}(U_{\rm
eff}(r,E^{\prime})+U_{\rm SO}(r,E^{\prime})\ {\bf L}\cdot{\bf
S})\\}\psi_{u}({\bf r})$
$\displaystyle=[(E^{\prime}-V_{c}(r))^{2}-m^{2}]\psi_{u}({\bf r})\,,$ (26)

where ${\bf S}$ and ${\bf L}$ are the total spin and angular momentum of the
nucleus, respectively. We can identify $U_{\rm eff}$ and $U_{\rm SO}$ as
effective Schroedinger-equation central and spin-orbit potentials that can be
constructed from $U_{s}(r,E^{\prime})$ and $U_{0}(r,E^{\prime})$. Note that
the Schroedinger-equation central potential also includes the Darwin term
accounting for relativistic corrections, and its effect is more pronounced in
the nuclear interior Arnold  _et al._ (1981). The spin-orbit term is
significantly smaller than the central one, and for this reason it has been
neglected in the present work. In the remainder of this paper, we will denote
the potential obtained from Ref. Cooper  _et al._ (2009) retaining only the
central contribution as the Schroedinger potential.

There are two ways that the potential plays a role within the cascade
algorithm. Firstly, the potential modifies the hard interactions that occur
between nucleons, often referred to as in-medium modifications in the
literature. In this work, we only consider the non-relativistic in-medium
corrections as implemented in Ref. Pandharipande and Pieper (1992). To account
for in-medium corrections due to the nuclear potential, we modify the
differential cross section using

$\displaystyle\frac{d\sigma^{\prime}}{d\Omega}$
$\displaystyle=\frac{\lvert\mathbf{p}^{\prime}_{1}-\mathbf{p}^{\prime}_{2}\rvert}{m}\left\lvert\frac{\mathbf{p}^{\prime}_{1}}{m^{*}(p^{\prime}_{1},\rho)}-\frac{\mathbf{p}^{\prime}_{2}}{m^{*}(p^{\prime}_{2},\rho)}\right\rvert^{-1}$
$\displaystyle\times\frac{m^{*}\left(\sqrt{(p_{3}^{\prime 2}+p_{4}^{\prime
2})/2},\rho\right)}{m}\frac{d\sigma}{d\Omega},$ (27)

where $p^{\prime}_{1},p^{\prime}_{2}$ are the momenta of the incoming
propagating nucleons, and $p^{\prime}_{3},p^{\prime}_{4}$ are the momenta of
the outgoing propagating nucleons. The effective nucleon mass $m^{*}$ is given
as

$m^{*}(p^{\prime},\rho)=p^{\prime}\left(\frac{p^{\prime}}{m}+\frac{dU(p^{\prime},\rho)}{dp^{\prime}}\right)^{-1}.$
(28)

This in-medium correction approximates the in-medium matrix element to be the
same as the free matrix element, and that
$U(p^{\prime}_{1},\rho)+U(p^{\prime}_{2},\rho)\approx U(\sqrt{(p_{1}^{\prime
2}+p_{2}^{\prime 2})/2},\rho)$. We leave the expansion to the relativistic
case to a future work. Note that we assume the potential to remain the same
regardless of INC dynamics. While this is certainly an approximation that will
fail when the nucleus suffers a “hard” breakdown, it should be reasonable when
the number of exiting nucleons is much lower than the number of a nucleons in
the nucleus.

We also consider the long-distance effect of a background potential on the
nucleon as it propagates through the nucleus. We simulate a particle
propagating by classical Hamiltonian evolution of the system. The equations of
motion can be written as

$\displaystyle\frac{{\rm d}\mathbf{p}}{{\rm d}t}$
$\displaystyle=-\frac{\partial
H}{\partial\mathbf{q}}=\left(\frac{(U_{s}+m)}{\sqrt{\mathbf{p}^{2}+(U_{s}+m)^{2}}}\frac{\partial
U_{s}}{\partial|\mathbf{q}|}+\frac{\partial
U_{0}}{\partial|\mathbf{q}|}\right)\hat{\mathbf{q}},$ $\displaystyle\frac{{\rm
d}\mathbf{q}}{{\rm d}t}$ $\displaystyle=\frac{\partial
H}{\partial\mathbf{p}}=\left(\frac{(U_{s}+m)\frac{\partial
U_{s}}{\partial|\mathbf{p}|}+|\mathbf{p}|}{\sqrt{\mathbf{p}^{2}+(U_{s}+m)^{2}}}+\frac{\partial
U_{0}}{\partial|\mathbf{p}|}\right)\hat{\mathbf{p}}\,.$ (29)

The equations above are clearly a set of coupled differential equations. In
order to maintain conservation of energy, a symplectic integrator is used for
the evolution. Since these differential equations are coupled, traditional
symplectic integrators will not work. It was shown in Ref. Tao (2016) that it
is possible to use symplectic integrators by working with an augmented
Hamiltonian in an extended phase space. App. B provides technical details.

##  V Comparison with exclusive observables

We now proceed to the analysis of exclusive observables in electron-carbon
scattering. Exclusive quantities are particularly relevant for neutrino
experiments, especially those based on the liquid argon time projection
chamber (LArTPC) technology such as the SBN detectors Antonello  _et al._
(2015) or the future DUNE experiment Acciarri  _et al._ (2015). LArTPCs are
able to identify and reconstruct tracks of all charged particles in a neutrino
scattering event in exceptional detail. This capability allows these detectors
to reject backgrounds and optimize searches more efficiently. If we take as an
example the recent MicroBooNE search for single photons Abratenko  _et al._
(2021) as an explanation of the MiniBooNE low energy excess Aguilar-Arevalo
_et al._ (2021), we can appreciate the importance of exclusive quantities: the
one-photon-zero-proton sample has a background rate 7 times higher than the
one-photon-one-proton sample, and this can largely be attributed to the
inability to reconstruct the $\Delta\to N\gamma$ invariant mass in the absence
of a proton track. Several other examples can be made, but the point is that
describing correctly exclusive observables will be crucial in current and
future neutrino experiments.

The CLAS and e4v collaborations have recently reported a study of energy
reconstruction in electron-nucleus scattering data, using methods employed in
neutrino experiments Khachatryan  _et al._ (2021). The collaborations analyzed
electron scattering data taken with CLAS at JLab for three different beam
energies: $1.159$, $2.257$, and $4.453~{}\text{GeV}$. The detection thresholds
for hadrons were similar to thresholds at current and future neutrino
experiments. The analysis focused on the reconstruction of several exclusive
and differential quantities, such as incoming electron energy reconstruction
for 0$\pi$ events, calorimetric reconstructed energies for 1p0$\pi$ events,
transverse variables, proton multiplicity, and so on.

In what follows, we describe the comparison between our generator and CLAS
data, for all available observables, focusing on quasielastic electron-carbon
scattering. The CLAS/e4v collaborations have reweighted their data by a factor
$Q^{4}/{\rm GeV}^{4}$, where $Q^{2}=-q^{2}$ is the four-momentum transfer
which can be obtained with final and initial electron kinematics. This was
done to have a better comparison with neutrino events: at these energies,
while electron-nucleus scattering is dominated by photon exchange, neutrino-
nucleus scattering can be very well approximated by a four-fermion
interaction. Here, we do the same in an event by event basis.

We adopt the same CLAS acceptances and mimic the energy resolution as
described in Ref. Khachatryan  _et al._ (2021), after corrections for
undetected particles. The electron and proton energies are smeared by
$1.5(0.5)\%$ and $3(1)\%$ for the 1.159 (2.257 and 4.453) GeV beams,
respectively. Protons were detected with momentum $p_{p}>300~{}\text{MeV}$ and
angle with respect to the beam direction $12^{\circ}<\theta_{p}$. Electrons
were detected with energy $E_{e}>0.4,\,0.55,$ and $1.1~{}\text{GeV}$ for
$E_{\rm beam}=1.159,\,2.257$ and $4.453~{}\text{GeV}$, as well as angles with
respect to the beam direction

$\theta_{e}^{i}>\theta_{0}^{i}+\frac{\theta_{1}^{i}}{p_{e}~{}[{\rm GeV}]},$
(30)

where $p_{e}$ is the electron momentum, $i=1,2,3$ refers to the three beam
energies in increasing order,
$\theta_{0}^{i}=17^{\circ},\,16^{\circ},\,13.5^{\circ}$ and
$\theta_{1}^{i}=7^{\circ},\,10.5^{\circ},\,15^{\circ}$. Since we do not
simulate production and propagation of pions, we do not list their acceptances
here.

We start with the double differential cross section
$d^{2}\sigma/d\Omega_{e}dE_{e}$, where $\Omega_{e}$ is the solid angle and
$E_{e}$ is the outgoing electron energy, as a function of the energy transfer
$\omega\equiv E_{e}^{\rm in}-E_{e}$, for fixed outgoing electron angle of
$37.5^{\circ}$ with respect to the beam axis and beam energy $E_{e}^{\rm
in}=1.159$ GeV, see Fig. 3. Hereafter we present ACHILLES results for several
different variations on the implementation of the INC, namely, the nucleon-
nucleon interaction model (Cylinder vs. Gaussian, see Ref. Isaacson  _et al._
(2021a)) and the real part of the nuclear potential (WFF, Schroedinger or
none). Here “none,” is used as a baseline prediction for the model described
in Ref. Isaacson  _et al._ (2021a). Different treatments will be color coded
and indicated by an inset in all figures. The spread among the lines can be
interpreted as one of the theoretical uncertainties on the lepton-nucleus
interaction modeling. Inclusive observables, such as those displayed in Fig. 3
are not affected by the semi-classical intranuclear cascades. Therefore, the
different lines lie on top of each other.

At low energy transfer, quasielastic scattering dominates the cross section.
In this region, particularly for $0.1<\omega<0.4$ GeV, our generator describes
the data fairly well, except for the small-energy region where the theory
underestimates the data. The agreement with data would be improved by the
interference effects neglected in intranuclear cascades (as discussed at the
end of Sec. III.2), yielding an enhancement of the strength at low $\omega$.
However, a naive combination of our intranuclear cascade and a folding
function would result in a double-counting for exclusive observables, as
discussed above. For this reason it has not been included in our calculation.

The missing strength in the quasielastic region and towards higher energy
transfers is largely ascribed to meson exchange, resonance production and deep
inelastic scattering contributions currently neglected in our analysis Rocco
_et al._ (2019). As one goes beyond this region towards higher energy
transfers, the quasielastic contribution shrinks and one expects other
components of the cross section to be more relevant, which explains the
discrepancy between the data and our generator. Overall, this level of
agreement is an encouraging result.

Figure 3: Comparison of the inclusive cross-section for an electron beam of
1159 MeV with outgoing angle of 37.5∘. Data is taken from Ref. Khachatryan
_et al._ (2021).

Another comparison we make is on the lepton energy reconstruction assuming
quasielastic scattering. The quasielastic energy reconstruction is done based
off the methodology used by water Cherenkov detectors, such as MiniBooNE and
T2K. In this case, only charged leptons and pions are measured. Assuming that
the neutrino scatters quasielastically from a stationary nucleon within a
nucleus, its incoming energy can be reconstructed as

$E_{\rm
QE}=\frac{2m_{N}\epsilon+2m_{N}E_{\ell}-m_{\ell}^{2}}{2\left(m_{N}-E_{\ell}+p_{\ell}\cos\theta_{\ell}\right)},$
(31)

where $m_{N}$ is the mass of the nucleon, $\epsilon$ is the average nucleon
separation energy (we use 21 MeV for carbon), $E_{\ell}(p_{\ell})$ is the
energy (momentum) of the outgoing lepton, and $\theta_{\ell}$ is the angle of
the outgoing lepton with respect to the beam axis. The different scheme
choices discussed in this paper are compared to the measured $E_{\rm QE}$
distribution for a 1.159 GeV electron beam on carbon from the CLAS data
Khachatryan  _et al._ (2021) in Fig. 4.

Here the peak around the beam energy is dominated by the quasielastic
contribution, while the tail towards lower values of $E_{\rm QE}$ is dominated
by meson exchange currents and resonance production. Therefore, we only expect
our results to approximately reproduce the peak, which is what is shown. The
agreement with the data for larger values of $E_{\rm QE}$ is likely to be
improved by the interference effects neglected by intranuclear cascades.
However, a more detailed analysis of the discrepancy will be carried out in
the future when meson exchange currents are included in ACHILLES. Analogously
to Fig. 3, this distribution has no information about the outgoing protons
contained within it. Therefore, we expect that the prediction should be
insensitive to the cascade parameters, as can be seen in the small spread of
the colored lines.

Figure 4: Comparison of the quasielastic energy reconstructed for an electron
beam of 1159 MeV. Data is taken from Ref. Khachatryan  _et al._ (2021). The
definition of $E_{QE}$ can be found in Eq. 31. The red dashed vertical line
marks the true beam energy.

In liquid argon time projection chamber experiments, such as MicroBooNE and
DUNE, the ionization energy is used as a mean to reconstruct the incoming
neutrino energy. In this case, the calorimetric energy is defined as

$E_{\rm cal}=\sum_{i}(E_{i}+\epsilon_{i}),$ (32)

where $E_{i}$ is the energy of the lepton or pions or the kinetic energy of
the protons, and $\epsilon$ is the average nucleon separation energy. In the
CLAS data, $E_{\rm cal}$ was calculated for events that contained exactly one
proton and zero pions Khachatryan  _et al._ (2021). The comparison between the
different schemes and the data are shown in Fig. 5, with beam energies of 1159
MeV in the top panel, 2257 MeV in the middle panel, and 4453 MeV in the bottom
panel. Since neutrons do not contribute to the calorimetric energy, we expect
this observable to be sensitive to the modeling of the intranuclear cascade.
The peak of these distributions correspond to the beam energy and is dominated
by the quasielastic contribution. The tail towards lower energies is due to
the intranuclear cascade, a result of the proton interacting with other
nucleons as it escapes the nucleus, as well as non-quasielastic interactions,
which are not currently implemented in ACHILLES. Around the peak, the largest
difference in peak height due to distinct implementations of the INC is about
7%. To be conservative, we will quote the INC theory uncertainty as the
largest difference among all INC implementations.

Figure 5: Comparison of the calorimetric energy reconstructed for an electron
beam of 1159 MeV (top), 2257 MeV (middle), and 4453 MeV (bottom). Data is
taken from Ref. Khachatryan  _et al._ (2021). The definition of $E_{cal}$ can
be found in Eq. (32).

To further study the accuracy of event simulation, the CLAS and e4v
collaborations studied three different transverse momentum related
observables. The first observable is the transverse momentum defined as

$\mathbf{p}^{T}=\mathbf{p}^{T}_{e}+\mathbf{p}^{T}_{p},$ (33)

where $\mathbf{p}^{T}_{e,p}$ is the transverse vector momentum with respect to
the beam axis for the electron and proton, respectively. Note that a lepton
scattering off protons at rest would only lead to
$p^{T}\equiv|\mathbf{p}^{T}|=0$. Fermi motion of nucleons in the nucleus will
lead to a distribution of $p^{T}$ around 100-200 MeV, while the intranuclear
cascade can introduce a long tail towards large values of $p^{T}$. This
observable is compared in Fig. 6 to the e4v data. The cascade tends to broaden
the spectrum in the quasielastic region, increasing the maximum value
observed. The spread due to different INC implementations is about 6%. The
$p^{T}>0.2$ GeV region has significant contributions from non-quasielastic
processes.

Figure 6: Comparison of the perpendicular momentum for an electron beam of
2257 MeV. Data is taken from Ref. Khachatryan  _et al._ (2021). The definition
of $p^{T}$ can be found in Eq. 33. The vertical dashed red lines denote the
location of the cuts on $p^{T}$ at 200 MeV and 400 MeV.

To isolate contributions from different nuclear processes, a cut is applied in
the $p^{T}$ variable before constructing the $E_{cal}$ distributions. The
results are shown in Fig. 7 for a cut of $p^{T}<200~{}\text{MeV}$ (top panel),
$200~{}\text{MeV}<p^{T}<400~{}\text{MeV}$ (middle panel), and
$p^{T}>400~{}\text{MeV}$ (bottom panel). Again, intranuclear cascades affect
the low $E_{cal}$ tail significantly, together with non-quasielastic
interactions. This is most evident in the $p^{T}>400$ MeV plot, in which both
effects are expected to have large impact. We find a theory uncertainty
associated to the INC implementation of 5% to 6% near the peak of all
distributions.

Figure 7: Comparison of the calorimetric energy reconstructed for an electron
beam of 2257 MeV, with a cut on the perpendicular momentum of $p^{T}<200$ MeV
(top), 200 MeV $<p^{T}<400$ MeV (middle), and $p^{T}>400$ MeV (bottom). Data
is taken from Ref. Khachatryan  _et al._ (2021).

The other two observables we use to validate ACHILLES are111Note there are
missing minus signs in Eqs. (7) and (8) of Ref. Khachatryan  _et al._ (2021)

$\displaystyle\delta\alpha_{T}$
$\displaystyle=\arccos\frac{-\mathbf{p}^{T}_{e}\cdot\mathbf{p}^{T}}{p^{T}_{e}p^{T}},$
(34) $\displaystyle\delta\phi_{T}$
$\displaystyle=\arccos\frac{-\mathbf{p}^{T}_{e}\cdot\mathbf{p}^{T}_{p}}{p^{T}_{e}p^{T}_{p}}.$
(35)

Note that $\mathbf{p}^{T}=-\mathbf{q}^{T}$, so $\delta\alpha_{T}$ is the angle
between the overall transverse momentum and the transverse momentum transfer.
Our results for $\delta\alpha_{T}$ are shown in Fig. 8. In the limit of no
final state interactions, $\mathbf{p}^{T}$ is simply the initial proton
transverse momentum. Since the initial proton momentum is isotropic,
$\delta\alpha_{T}$ should also be isotropic in this limit. The increase in the
high-angle region of the $\delta\alpha_{T}$ distribution can be attributed to
intranuclear cascades and non-quasielastic interactions. We find an INC theory
uncertainty in $\delta\alpha_{T}$ of about 10%.

On the other hand, $\delta\phi_{T}$ measures the opening angle between the
proton and the transverse momentum transfer. We present the comparison to data
in Fig. 9. In the absence of both final state interactions and initial proton
momentum, we have $\mathbf{p}^{T}_{p}=-\mathbf{p}^{T}_{e}$ and thus
$\delta\phi_{T}$ is a delta function at zero. In the presence of Fermi
momentum, the struck proton has a nonzero momentum, $k_{p}\neq 0$, which
smears the $\delta\phi_{T}$ distribution around zero by
$\mathcal{O}(k^{T}_{p}/p^{T}_{e})$. Final state interactions help to smear out
the distribution to larger opening angles, partially explaining the high
$\delta\phi_{T}$ tails in Fig. 9. The INC uncertainty is found to be about 5%.

Finally, notice that observables for which final-state interactions play an
important role offer the greatest sensitivity to the implementation of the INC
model. This sensitivity is visible in the spread in the color histograms and
is particularly evident, for example, in the $E_{cal}$ distribution for higher
$p_{T}$ (see the bottom panel of Fig. 7) and in the high angle region of the
$\alpha_{T}$ distribution (see top and middle panels in Fig. 8).

Figure 8: Comparison of $\delta\alpha_{T}$ for an electron beam of 1159 MeV
(top), 2257 MeV (middle), 4453 MeV (bottom). Data is taken from Ref.
Khachatryan  _et al._ (2021). The definition of $\delta\alpha_{T}$ can be
found in Eq. 34.

Figure 9: Comparison of $\delta\phi_{T}$ for an electron beam of 1159 MeV
(top), 2257 MeV (middle), 4453 MeV (bottom). Data is taken from Ref.
Khachatryan  _et al._ (2021). The definition of $\delta\phi_{T}$ can be found
in Eq. 35.

##  VI Other observables

In this section we propose additional key observables that could be measured
in current and future electron-nucleus scattering experiments, such as CLAS12
Burkert  _et al._ (2020) or LDMX Åkesson  _et al._ (2018b). The goal is to
encourage the experimental collaborations to present such observables that
will ultimately serve to validate lepton-nucleus interaction models.

Let us start with an exclusive differential observable that is highly
sensitive to final state interactions: the _proton multiplicity energy
spectrum_. As we currently do not have pions propagating in our intranuclear
cascade modeling in ACHILLES, we focus on $n$p0$\pi$ events. Taking the 2.257
GeV electron beam as an example, for every event, we count the number of
protons that pass experimental cuts (see Sec. V). Then we take all leading-
energy protons in events with at least one proton and build their energy
spectrum. We repeat the procedure for all second- and third-leading protons,
in events with at least two or three protons, respectively.

The results of this procedure are the proton energy spectra shown in Fig. 10,
from the leading proton in the upper panel to the third leading proton in the
lower panel. We would expect this distribution to be highly sensitive to the
intranuclear cascade model. INCs may raise the proton multiplicity,
contributing to the spectra of second and third protons, and tend to
distribute the energy among all outgoing protons, shifting the leading proton
spectrum towards lower energies. This is indeed observed when comparing the
INC uncertainties in peak regions in the three panels of Fig. 10, which are
approximately 3%, 12%, and 15%, from top to bottom. We also expect that other
cross section channels will significantly contribute to this observable in a
nontrivial way. For example, while DIS occurs for higher momentum transfer,
hadronization followed by final state interactions may lead to a large
multiplicity of low energy protons. The proton multiplicity energy spectra for
all interaction channels will be a subject of a future publication.

Figure 10: Energy spectra of the $n$-th most energetic proton, from the most
(upper panel) to the third most energetic protons (lower panel) for an
electron beam with an energy of 2.257 GeV.

In the same vein, we also propose the _proton multiplicity angular spectra_.
We take again the 2.257 GeV electron beam as an example. We plot the angle of
the leading, second and third protons with respect to the beam axis in Fig.
11, from the leading proton in the upper panel to the third leading proton in
the lower panel. Note that we have decided to order the protons according to
their energies. Our main motivation lies on the fact that higher energy
protons are more relevant to the reconstruction of neutrino energies, and
therefore a correct description of the leading protons is more relevant than
the subleading ones. Again, we expect intranuclear cascade models to play a
crucial role here, as well as the other interactions channels, which will be
the studied in a future publication.

Figure 11: Angular spectra with respect to the beam axis of the $n$-th most
energetic proton, from the most (upper panel) to the third most energetic
protons (lower panel) for an electron beam with an energy of 2.257 GeV.

Another interesting observable is the angle between the sum of the momenta of
all visible outgoing particles with respect to the beam axis. The only
particle we take to be invisible here are neutrons. We apply the 1p0$\pi$
selection cuts from the CLAS/e4v analysis, see Sec. V. This angle would be
zero in the case of an electron scattering on a free proton at rest. This
observable is motivated the physics of atmospheric neutrinos. In this sample,
the incoming neutrino direction needs to be reconstructed in order to estimate
the neutrino path through the Earth and in the oscillation probabilities. A
measurement of atmospheric neutrinos in the $0.1-1$ GeV scale at the DUNE
experiment could provide nontrivial information on the $CP$ violation phase
Kelly  _et al._ (2019), and could also be used to perform a tomography study
of the Earth, contributing to our understanding of the chemical composition of
its core Kelly  _et al._ (2021).

To be more precise, we define the reconstructed beam angle in electron-nucleus
scattering as

$\cos\theta_{\rm rec}\equiv\frac{\hat{\mathbf{k}}_{e}\cdot\mathbf{p}_{\rm
out}}{|\mathbf{p}_{\rm out}|},$ (36)

where $\mathbf{p}_{\rm out}$ is the sum of all momenta of visible outgoing
particles and $\hat{\mathbf{k}}_{e}$ is a unit vector in the beam direction.
The reconstructed angle $\theta_{\rm rec}$ can deviate from zero for several
reasons: Fermi motion, as it adds momentum to the incoming proton that is not
accounted for in Eq. (36); intranuclear cascade, as a proton may scatter off a
neutron which in turn may be invisible to most detectors of interest; and
nuclear potential, which may deflect the outgoing proton. Our results are
found in Fig. 12. The spread around $\theta_{\rm rec}=0$ can be attributed to
Fermi motion and is of order $\Delta\theta_{\rm rec}\sim p_{F}/k_{e}$, where
$p_{F}$ is the Fermi momentum. The tail at large angles is due to final state
interactions, as protons may scatter off neutrons and deflect significantly.
Non-quasielastic interactions should further populate the high $\theta_{\rm
rec}$ region.

Figure 12: Reconstructed angle of incoming electron with respect to the beam
axis for an electron beam with an energy of 2.257 GeV.

##  VII Conclusions

We have presented a newly developed lepton-nucleus event generator, ACHILLES.
Our generator factorizes the primary interaction vertex from the propagation
of hadrons througout the nucleus, allowing for a great deal of modularity,
which is one of the pillars of ACHILLES. Due to this modularity, ACHILLES can
be used for generating either electron-nucleus or neutrino-nucleus scattering
events, and the implementation of numerous scenarios for physics beyond the
Standard Model is straightforward.

We have validated quasielastic scattering against high quality, inclusive and
exclusive, electron-carbon data, including the recent CLAS/e4v reanalysis of
existing data. We find good agreement between data and simulation. A complete
estimate of the theoretical uncertainty associated with the nuclear model and
the current operator adopted in the description of the primary interaction
vertex is highly non-trivial and has not been included in this work. A
promising avenue to quantify model dependence involves testing different
nuclear many-body methods, possibly including different nuclear currents, form
factors, and Hamiltonians as inputs. A study along these lines has been
carried out in Ref. Andreoli  _et al._ (2022) where the inclusive differential
cross sections for electron scattering on 3He and 3H have been evaluated using
different many-body approaches based on the same description of nuclear
dynamics. Inputs from LQCD calculations, such as nucleon form factors and
elementary nucleon matrix elements, will be incorporated as they become
available.

By varying model assumptions of the intranuclear cascade (namely, different
nucleon-nucleon interactions models and nuclear potentials), we have estimated
one component of the overall theory uncertainty budget in electron-nucleus
scattering. For observables that are sensitive to final-state interactions,
the theoretical model dependence associated with different intranuclear
cascade models is typically a 5-10% effect. Theory uncertainty estimates will
be crucial for a precision neutrino physics program, in particular for the
DUNE experiment.

We have also proposed novel observables that will allow for further validation
of lepton-nucleus scattering models. Although we have only analyzed electron-
carbon scattering data in the quasielastic region, our code is readily
extendable to generate neutrino-nucleus scattering events. Comparison against
neutrino scattering data, as well as the inclusion of other primary
interaction modes, such as resonant scattering, meson exchange current and
deep inelastic scattering, will be subjects of future publications.

##  VIII Acknowledgments

We thank Or Hen and Afroditi Papadopoulou for useful discussions. We thank Adi
Ashkenazi, Or Hen, Stefan Höche, Shirley Li, Kendall Mahn, Afroditi
Papadopoulou, Luke Pickering, and Larry Weinstein for comments on the
manuscript. Fermilab is operated by the Fermi Research Alliance, LLC under
contract No. DE-AC02-07CH11359 with the United States Department of Energy.
The present research is supported by the U.S. Department of Energy, Office of
Science, Office of Nuclear Physics, under contracts DE-AC02-06CH11357 (A.L).
A.L acknowledges funding from the INFN grant INNN3, and from the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No 824093. This material is based upon work supported by the U.S. Department
of Energy, Office of Science, Office of Nuclear Physics under grant Contract
Numbers DE-SC0011090 and DE-SC0021006 (W.J.). This project has received
support from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 860881-HIDDeN.

##  Appendix A ACHILLES Technical Details

The calculation of the primary interaction within the ACHILLES generator is
separated into a leptonic and a hadronic current as described in Isaacson  _et
al._ (2021b). The use of currents reduces the bookkeeping required to properly
handle the interference between different gauge bosons contributing in the
primary interaction. This is important when dealing with BSM scenarios in
which the dominant contribution to the total cross section may arise from the
interference with the Standard Model. The setup of initial and final state
particles, along with the configuration of all other options within the
generator is controlled with a set of YAML files. Almost all parameters can be
controlled through the YAML input card without the need to recompile the code.

While the core of ACHILLES is written in modern C++ for high performance, a
general purpose fortran90 wrapper is provided to interface ACHILLES to
available nuclear models. The wrapper consists of three components. Firstly,
the ACHILLES code provides an interface to the physical constants and other
useful common utilities, such as the handling of particle information
(particle id, mass, four-momentum, position, status code, etc.). This helps to
ensure consistency of physical constants and particles throughout the
calculation. Secondly, the ACHILLES code expects the nuclear model to define
two functions that define the interface between the C++ and fortran90. The
first function handles the initialization of the nuclear model, which gets
passed as an argument, a filename and length to be handled by the fotran90
code. The second function is expected to perform the calculation of the
nuclear current, as discussed in the previous paragraph. This function is
passed as input information about the four momentum of all the nucleons and
the gauge boson for a given event, and expects to be returned the nuclear
current. Finally, the nuclear model needs to be registered with the ACHILLES
code to provide a means to enable simulations of the model via the input card.
This wrapper was used to include the original, extensively validated spectral
function codes written in fortran90 to be used for the nuclear initial state
into the ACHILLES generator.

The sampling of the phase space is performed with the efficient multi-channel
Kleiss and Pittau (1994) and recursive phase space Byckling and Kajantie
(1969) discussed in Isaacson  _et al._ (2021b) with importance sampling
handled by the VEGAS algorithm Lepage (1978, 1980). This makes the code
readily extendable to other reaction mechanisms and higher dimensional phase
spaces. Based on the experience of the LHC community, we do not expect any
dramatic decrease in computational speed as the multiplicity of the final
state increases beyond the extra time involved in evaluating the matrix
element and generating the additional momentum. Finally, the major benefit in
using these sampling techniques is in the increased unweighting efficiency
during the event generation process.

Further details on the input card and API details on the interface will be
expanded upon in a manual to be released in the future.

##  Appendix B Symplectic Integrator

As detailed in Ref. Tao (2016), in order to develop an explicit symplectic
integrator for non-separable Hamilitonians can be achieved by using an
augmented Hamiltonian defined as follows:

$\overline{H}(q,p,x,y)\equiv H_{A}(q,y)+H_{B}(x,p)+\omega H_{C}(q,p,x,y),$
(37)

where $H_{A}(q,y)$ and $H_{B}(x,p)$ are copies of the original Hamiltonian,
$H_{C}(q,p,x,y)=|q-x|^{2}/2+|p-y|^{2}/2$ acts as a harmonic oscillator keeping
the two solutions close in phase space, and $\omega$ is a tunable parameter to
control the strength of the coupling. Explicit flows can then be defined as:

$\displaystyle\phi_{H_{A}}^{\delta}$ $\displaystyle:\begin{bmatrix}q\\\ p\\\
x\\\ y\end{bmatrix}\to\begin{bmatrix}q\\\ p-\delta\partial_{q}H(q,y)\\\
x+\delta\partial_{y}H(q,y)\\\ y\end{bmatrix},$ (38)
$\displaystyle\phi_{H_{B}}^{\delta}$ $\displaystyle:\begin{bmatrix}q\\\ p\\\
x\\\ y\end{bmatrix}\to\begin{bmatrix}q+\delta\partial_{p}H(x,p)\\\ p\\\ x\\\
y-\delta\partial_{x}H(x,p)\end{bmatrix},$ (39)
$\displaystyle\phi_{H_{C}}^{\delta}$ $\displaystyle:\begin{bmatrix}q\\\ p\\\
x\\\ y\end{bmatrix}\to\frac{1}{2}\begin{bmatrix}\begin{pmatrix}q+x\\\
p+y\end{pmatrix}+R(\delta)\begin{pmatrix}q-x\\\ p-y\end{pmatrix}\\\
\begin{pmatrix}q+x\\\ p+y\end{pmatrix}-R(\delta)\begin{pmatrix}q-x\\\
p-y\end{pmatrix}\end{bmatrix},$ (40)

where $\delta$ is the time step for the evolution, and

$R(\delta)\equiv\begin{bmatrix}\cos(2\omega\delta)&\sin(2\omega\delta)\\\
-\sin(2\omega\delta)&\cos(2\omega\delta)\end{bmatrix}.$ (41)

A second order method can be created from these three Hamiltonians through the
use of the symmetric Strang splitting method Strang (1968). This leads to a
single evolution step over a time step $\delta$ as:

$\phi_{2}^{\delta}=\phi_{H_{A}}^{\delta/2}\circ\phi_{H_{B}}^{\delta/2}\circ\phi_{\omega
H_{C}}^{\delta}\circ\phi_{H_{B}}^{\delta/2}\circ\phi_{H_{A}}^{\delta/2}.$ (42)

A $l$th order integration method can be obtained through the triple jump
method Yoshida (1990); McLachlan and Quispel (2002):

$\displaystyle\phi_{l}^{\delta}$
$\displaystyle=\phi_{l-2}^{\gamma_{l}\delta}\circ\phi_{l-2}^{(l-2)\gamma_{l}\delta}\circ\phi_{l-2}^{\gamma\delta},$
(43) $\displaystyle\text{where}\ \gamma_{l}$
$\displaystyle=\frac{1}{2-2^{1/(l+1)}},$

which will also be symplectic if $\phi_{l-2}$ is symplectic. In this work, we
tune the values of $\omega$ and $\delta$ such that the results are stable as a
second order integrator.

Figure 13: Demonstration of the stability of the symplectic integrator. The
left panel shows the trajectory of a nucleon with a momentum of 250 MeV
perpendicular to the radius starting at a radius of $r=1$ fm in the Wiringa
potential (blue) and the Cooper potential (red) over 100,000 time steps. The
right panel shows the deviation from the starting energy as a function of the
time step.

Figure 13 demonstrates the stability of the symplectic integrator for a
nucleon with a momentum of 250 MeV perpendicular to the radius starting at a
radius of $r=1$ fm in the non-relativisitc Wiringa potential (blue) and the
relativistic Cooper potential (red). The simulation is run for 100,000 time
steps, and the maximum energy deviation is of the order of $10^{-4}$. The
deviation is periodic in nature, which is a common feature for symplectic
integrators.

## References

  * Itow  _et al._ (2001) Y. Itow _et al._ (T2K), in  _3rd Workshop on Neutrino Oscillations and Their Origin (NOON 2001)_ (2001) arXiv:hep-ex/0106019 .
  * Ayres  _et al._ (2004) D. S. Ayres _et al._ (NOvA), (2004), arXiv:hep-ex/0503053 .
  * Acciarri  _et al._ (2015) R. Acciarri _et al._ (DUNE), (2015), arXiv:1512.06148 [physics.ins-det] .
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Abbott  _et al._ (1998) D. Abbott _et al._ , Phys. Rev. Lett. 80, 5072 (1998).
  * Dutta  _et al._ (2000) D. Dutta _et al._ , Phys. Rev. C 61, 061602 (2000).
  * Mecking  _et al._ (2003) B. A. Mecking _et al._ (CLAS), Nucl. Instrum. Meth. A 503, 513 (2003).
  * Dutta  _et al._ (2003) D. Dutta _et al._ (JLab E91013), Phys. Rev. C 68, 064603 (2003), arXiv:nucl-ex/0303011 .
  * Rohe  _et al._ (2005) D. Rohe _et al._ (E97-006), Phys. Rev. C 72, 054602 (2005), arXiv:nucl-ex/0506007 .
  * Burkert  _et al._ (2020) V. D. Burkert _et al._ , Nucl. Instrum. Meth. A 959, 163419 (2020).
  * Battaglieri  _et al._ (2017) M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter_ (2017) arXiv:1707.04591 [hep-ph] .
  * Kahn  _et al._ (2018) Y. Kahn, G. Krnjaic, N. Tran, and A. Whitbeck, JHEP 09, 153 (2018), arXiv:1804.03144 [hep-ph] .
  * Åkesson  _et al._ (2018a) T. Åkesson _et al._ (LDMX), (2018a), arXiv:1808.05219 [hep-ex] .
  * Battaglieri  _et al._ (2019) M. Battaglieri _et al._ (BDX), (2019), arXiv:1910.03532 [physics.ins-det] .
  * Agrawal  _et al._ (2021) P. Agrawal _et al._ , (2021), arXiv:2102.12143 [hep-ph] .
  * Battaglieri  _et al._ (2021) M. Battaglieri _et al._ , (2021), arXiv:2105.04540 [hep-ex] .
  * Serber (1947) R. Serber, Phys. Rev. 72, 1114 (1947).
  * Metropolis  _et al._ (1958) N. Metropolis, R. Bivins, M. Storm, J. M. Miller, G. Friedlander, and A. Turkevich, Phys. Rev. 110, 204 (1958).
  * Bertini (1963) H. W. Bertini, Phys. Rev. 131, 1801 (1963).
  * Cugnon (1980) J. Cugnon, Phys. Rev. C 22, 1885 (1980).
  * Bertsch  _et al._ (1984) G. F. Bertsch, H. Kruse, and S. D. Gupta, Phys. Rev. C 29, 673 (1984), [Erratum: Phys.Rev.C 33, 1107–1108 (1986)].
  * Stoecker and Greiner (1986) H. Stoecker and W. Greiner, Phys. Rept. 137, 277 (1986).
  * Bauer  _et al._ (1986) W. Bauer, G. F. Bertsch, W. Cassing, and U. Mosel, Phys. Rev. C 34, 2127 (1986).
  * Bertsch and Das Gupta (1988) G. F. Bertsch and S. Das Gupta, Phys. Rept. 160, 189 (1988).
  * Botermans and Malfliet (1990) W. Botermans and R. Malfliet, Phys. Rept. 198, 115 (1990).
  * Cassing  _et al._ (1990) W. Cassing, V. Metag, U. Mosel, and K. Niita, Phys. Rept. 188, 363 (1990).
  * Danielewicz and Bertsch (1991) P. Danielewicz and G. F. Bertsch, Nucl. Phys. A 533, 712 (1991).
  * Teis  _et al._ (1997) S. Teis, W. Cassing, M. Effenberger, A. Hombach, U. Mosel, and G. Wolf, Z. Phys. A 356, 421 (1997), arXiv:nucl-th/9609009 .
  * Cugnon  _et al._ (1997) J. Cugnon, C. Volant, and S. Vuillier, Nucl. Phys. A 620, 475 (1997).
  * Golubeva  _et al._ (1998) Y. S. Golubeva, L. A. Kondratyuk, A. Bianconi, S. Boffi, and M. Radici, Phys. Rev. C 57, 2618 (1998), arXiv:nucl-th/9712040 .
  * Boudard  _et al._ (2002) A. Boudard, J. Cugnon, S. Leray, and C. Volant, Phys. Rev. C 66, 044615 (2002).
  * Hayato (2002) Y. Hayato, Nucl. Phys. B Proc. Suppl. 112, 171 (2002).
  * Casper (2002) D. Casper, Nucl. Phys. B Proc. Suppl. 112, 161 (2002), arXiv:hep-ph/0208030 .
  * Duarte (2007) H. Duarte, Phys. Rev. C 75, 024611 (2007).
  * Andreopoulos  _et al._ (2010) C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Iwamoto  _et al._ (2010) H. Iwamoto, M. Imamura, Y. Koba, Y. Fukui, G. Wakabayashi, Y. Uozumi, T. Kin, Y. Iwamoto, S. Hohara, and M. Nakano, Phys. Rev. C 82, 034604 (2010).
  * Buss  _et al._ (2012) O. Buss, T. Gaitanos, K. Gallmeister, H. van Hees, M. Kaskulov, O. Lalakulich, A. B. Larionov, T. Leitner, J. Weil, and U. Mosel, Phys. Rept. 512, 1 (2012), arXiv:1106.1344 [hep-ph] .
  * Golan  _et al._ (2012) T. Golan, C. Juszczak, and J. T. Sobczyk, Phys. Rev. C 86, 015505 (2012), arXiv:1202.4197 [nucl-th] .
  * Sawada  _et al._ (2012) Y. Sawada, Y. Uozumi, S. Nogamine, T. Yamada, Y. Iwamoto, T. Sato, and K. Niita, Nucl. Instrum. Meth. B 291, 38 (2012).
  * Uozumi  _et al._ (2012) Y. Uozumi, T. Yamada, S. Nogamine, and M. Nakano, Phys. Rev. C 86, 034610 (2012).
  * Battistoni  _et al._ (2013) G. Battistoni, F. Cerutti, A. Ferrari, J. Ranft, S. Roesler, and P. R. Sala, J. Phys. Conf. Ser. 408, 012051 (2013).
  * Battistoni  _et al._ (2015) G. Battistoni _et al._ , Annals Nucl. Energy 82, 10 (2015).
  * Mosel (2019) U. Mosel, J. Phys. G 46, 113001 (2019), arXiv:1904.11506 [hep-ex] .
  * Niewczas and Sobczyk (2019) K. Niewczas and J. T. Sobczyk, Phys. Rev. C 100, 015505 (2019), arXiv:1902.05618 [hep-ex] .
  * Isaacson  _et al._ (2021a) J. Isaacson, W. I. Jay, A. Lovato, P. A. N. Machado, and N. Rocco, Phys. Rev. C 103, 015502 (2021a), arXiv:2007.15570 [hep-ph] .
  * Dytman  _et al._ (2021) S. Dytman, Y. Hayato, R. Raboanary, J. T. Sobczyk, J. Tena Vidal, and N. Vololoniaina, Phys. Rev. D 104, 053006 (2021), arXiv:2103.07535 [hep-ph] .
  * Bhattacharya  _et al._ (2011) B. Bhattacharya, R. J. Hill, and G. Paz, Phys. Rev. D 84, 073006 (2011), arXiv:1108.0423 [hep-ph] .
  * Meyer  _et al._ (2016) A. S. Meyer, M. Betancourt, R. Gran, and R. J. Hill, Phys. Rev. D 93, 113015 (2016), arXiv:1603.03048 [hep-ph] .
  * Davoudi  _et al._ (2021) Z. Davoudi, W. Detmold, K. Orginos, A. Parreño, M. J. Savage, P. Shanahan, and M. L. Wagman, Phys. Rept. 900, 1 (2021), arXiv:2008.11160 [hep-lat] .
  * Borah  _et al._ (2020) K. Borah, R. J. Hill, G. Lee, and O. Tomalak, Phys. Rev. D 102, 074012 (2020), arXiv:2003.13640 [hep-ph] .
  * Meyer  _et al._ (2022) A. S. Meyer, A. Walker-Loud, and C. Wilkinson, (2022), arXiv:2201.01839 [hep-lat] .
  * Benhar  _et al._ (2008) O. Benhar, D. day, and I. Sick, Rev. Mod. Phys. 80, 189 (2008), arXiv:nucl-ex/0603029 .
  * Martini  _et al._ (2009) M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C 80, 065501 (2009), arXiv:0910.2622 [nucl-th] .
  * Martini  _et al._ (2010) M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C 81, 045502 (2010), arXiv:1002.4538 [hep-ph] .
  * Amaro  _et al._ (2011) J. E. Amaro, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, and C. F. Williamson, Phys. Lett. B 696, 151 (2011), arXiv:1010.1708 [nucl-th] .
  * Gran  _et al._ (2013) R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev. D 88, 113007 (2013), arXiv:1307.8105 [hep-ph] .
  * Benhar  _et al._ (2015) O. Benhar, A. Lovato, and N. Rocco, Phys. Rev. C 92, 024602 (2015), arXiv:1502.00887 [nucl-th] .
  * Lovato  _et al._ (2020) A. Lovato, J. Carlson, S. Gandolfi, N. Rocco, and R. Schiavilla, Phys. Rev. X 10, 031068 (2020), arXiv:2003.07710 [nucl-th] .
  * Sobczyk  _et al._ (2021) J. E. Sobczyk, B. Acharya, S. Bacca, and G. Hagen, Phys. Rev. Lett. 127, 072501 (2021), arXiv:2103.06786 [nucl-th] .
  * Roggero and Carlson (2019) A. Roggero and J. Carlson, Phys. Rev. C 100, 034610 (2019), arXiv:1804.01505 [quant-ph] .
  * Roggero  _et al._ (2020) A. Roggero, A. C. Y. Li, J. Carlson, R. Gupta, and G. N. Perdue, Phys. Rev. D 101, 074038 (2020), arXiv:1911.06368 [quant-ph] .
  * Kadanoff and Baym (1962) L. Kadanoff and G. Baym,  _Quantum Statistical Mechanics_ (Benjamin, New York, 1962).
  * Cugnon and Henrotte (2003) J. Cugnon and P. Henrotte, European Physical Journal A 16, 393 (2003).
  * Dore  _et al._ (2001) D. Dore _et al._ , Phys. Rev. C 63, 034612 (2001).
  * Buss  _et al._ (2008) O. Buss, T. Leitner, U. Mosel, and L. Alvarez-Ruso, in  _46th International Winter Meeting on Nuclear Physics_ (2008) arXiv:0809.2550 [nucl-th] .
  * Leitner  _et al._ (2009) T. Leitner, O. Buss, L. Alvarez-Ruso, and U. Mosel, Phys. Rev. C 79, 034601 (2009), arXiv:0812.0587 [nucl-th] .
  * Mosel and Gallmeister (2019) U. Mosel and K. Gallmeister, Phys. Rev. C 99, 064605 (2019), arXiv:1811.10637 [nucl-ex] .
  * Ankowski and Friedland (2020) A. M. Ankowski and A. Friedland, Phys. Rev. D 102, 053001 (2020), arXiv:2006.11944 [hep-ph] .
  * Khachatryan  _et al._ (2021) M. Khachatryan _et al._ (CLAS, e4v), Nature 599, 565 (2021).
  * Abi  _et al._ (2020a) B. Abi _et al._ (DUNE), JINST 15, T08008 (2020a), arXiv:2002.02967 [physics.ins-det] .
  * Abi  _et al._ (2020b) B. Abi _et al._ (DUNE), (2020b), arXiv:2002.03005 [hep-ex] .
  * Isaacson  _et al._ (2021b) J. Isaacson, S. Höche, D. L. Gutierrez, and N. Rocco, (2021b), arXiv:2110.15319 [hep-ph] .
  * Buckley  _et al._ (2021) A. Buckley, P. Ilten, D. Konstantinov, L. Lönnblad, J. Monk, W. Pokorski, T. Przedzinski, and A. Verbytskyi, Comput. Phys. Commun. 260, 107310 (2021), arXiv:1912.08005 [hep-ph] .
  * Buckley  _et al._ (2013) A. Buckley, J. Butterworth, D. Grellscheid, H. Hoeth, L. Lonnblad, J. Monk, H. Schulz, and F. Siegert, Comput. Phys. Commun. 184, 2803 (2013), arXiv:1003.0694 [hep-ph] .
  * Bierlich  _et al._ (2020) C. Bierlich _et al._ , SciPost Phys. 8, 026 (2020), arXiv:1912.05451 [hep-ph] .
  * Stowell  _et al._ (2017) P. Stowell _et al._ , JINST 12, P01016 (2017), arXiv:1612.07393 [hep-ex] .
  * Rocco  _et al._ (2019) N. Rocco, S. X. Nakamura, T. S. H. Lee, and A. Lovato, Phys. Rev. C 100, 045503 (2019), arXiv:1907.01093 [nucl-th] .
  * Ruso  _et al._ (2022) L. A. Ruso _et al._ , (2022), arXiv:2203.09030 [hep-ph] .
  * Carlson  _et al._ (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. 87, 1067 (2015), arXiv:1412.3081 [nucl-th] .
  * Höche (2015) S. Höche, in  _Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders_ (2015) pp. 235–295, arXiv:1411.4085 [hep-ph] .
  * Marcucci  _et al._ (2005) L. E. Marcucci, M. Viviani, R. Schiavilla, A. Kievsky, and S. Rosati, Phys. Rev. C 72, 014001 (2005), arXiv:nucl-th/0502048 .
  * Benhar and Pandharipande (1993) O. Benhar and V. R. Pandharipande, Phys. Rev. C 47, 2218 (1993).
  * Benhar  _et al._ (1994) O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A 579, 493 (1994).
  * Ivanov  _et al._ (2019) 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. Udías, Phys. Rev. C 99, 014610 (2019), arXiv:1812.09435 [nucl-th] .
  * Rocco and Barbieri (2018) N. Rocco and C. Barbieri, Phys. Rev. C 98, 025501 (2018), arXiv:1803.00825 [nucl-th] .
  * Barbieri  _et al._ (2019) C. Barbieri, N. Rocco, and V. Somà, Phys. Rev. C 100, 062501 (2019), arXiv:1907.01122 [nucl-th] .
  * Andreoli  _et al._ (2022) L. Andreoli, J. Carlson, A. Lovato, S. Pastore, N. Rocco, and R. B. Wiringa, Phys. Rev. C 105, 014002 (2022), arXiv:2108.10824 [nucl-th] .
  * Rocco (2020) N. Rocco, Front. in Phys. 8, 116 (2020).
  * Ruiz Simo  _et al._ (2017) I. Ruiz Simo, J. E. Amaro, M. B. Barbaro, A. De Pace, J. A. Caballero, and T. W. Donnelly, J. Phys. G 44, 065105 (2017), arXiv:1604.08423 [nucl-th] .
  * Hernandez  _et al._ (2007) E. Hernandez, J. Nieves, and M. Valverde, Phys. Rev. D 76, 033005 (2007), arXiv:hep-ph/0701149 .
  * Kamano  _et al._ (2013) H. Kamano, S. X. Nakamura, T. S. H. Lee, and T. Sato, Phys. Rev. C 88, 035209 (2013), arXiv:1305.4351 [nucl-th] .
  * Nakamura  _et al._ (2015) S. X. Nakamura, H. Kamano, and T. Sato, Phys. Rev. D 92, 074024 (2015), arXiv:1506.03403 [hep-ph] .
  * Kamano  _et al._ (2016) H. Kamano, S. X. Nakamura, T. S. H. Lee, and T. Sato, Phys. Rev. C 94, 015201 (2016), arXiv:1605.00363 [nucl-th] .
  * O’Connell  _et al._ (1987) J. S. O’Connell _et al._ , Phys. Rev. C 35, 1063 (1987).
  * Sealock  _et al._ (1989) R. M. Sealock _et al._ , Phys. Rev. Lett. 62, 1350 (1989).
  * Zeller (1973) D. Zeller,  _INVESTIGATION OF THE STRUCTURE OF THE C-12 NUCLEUS BY HIGH-ENERGY ELECTRON SCATTERING_ , Other thesis (1973).
  * Cooper  _et al._ (2009) E. D. Cooper, S. Hama, and B. C. Clark, Phys. Rev. C 80, 034605 (2009).
  * Benhar (2013) O. Benhar, Phys. Rev. C 87, 024606 (2013), arXiv:1301.3357 [nucl-th] .
  * González-Jiménez  _et al._ (2019) R. González-Jiménez, A. Nikolakopoulos, N. Jachowicz, and J. M. Udías, Phys. Rev. C 100, 045501 (2019), arXiv:1904.10696 [nucl-th] .
  * Ankowski  _et al._ (2015) A. M. Ankowski, O. Benhar, and M. Sakuda, Phys. Rev. D 91, 033005 (2015), arXiv:1404.5687 [nucl-th] .
  * Ershova  _et al._ (2022) A. Ershova _et al._ , (2022), arXiv:2202.10402 [hep-ph] .
  * Nikolakopoulos  _et al._ (2022) A. Nikolakopoulos, R. González-Jiménez, N. Jachowicz, K. Niewczas, F. Sánchez, and J. M. Udías, (2022), arXiv:2202.01689 [nucl-th] .
  * Wiringa (1988) R. B. Wiringa, Phys. Rev. C 38, 2967 (1988).
  * Wiringa  _et al._ (1988) R. B. Wiringa, V. Fiks, and A. Fabrocini, Phys. Rev. C 38, 1010 (1988).
  * Clark  _et al._ (2006) B. C. Clark, E. D. Cooper, and S. Hama, Phys. Rev. C 73, 024608 (2006).
  * Arnold  _et al._ (1981) L. G. Arnold, B. C. Clark, R. L. Mercer, and P. Schwandt, Phys. Rev. C 23, 1949 (1981).
  * Pandharipande and Pieper (1992) V. R. Pandharipande and S. C. Pieper, Phys. Rev. C 45, 791 (1992).
  * Tao (2016) M. Tao, Phys. Rev. E 94, 043303 (2016).
  * Antonello  _et al._ (2015) M. Antonello _et al._ (MicroBooNE, LAr1-ND, ICARUS-WA104), (2015), arXiv:1503.01520 [physics.ins-det] .
  * Abratenko  _et al._ (2021) P. Abratenko _et al._ (MicroBooNE), (2021), arXiv:2110.00409 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2021) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 103, 052002 (2021), arXiv:2006.16883 [hep-ex] .
  * Åkesson  _et al._ (2018b) T. Åkesson _et al._ (LDMX), (2018b), arXiv:1808.05219 [hep-ex] .
  * Kelly  _et al._ (2019) K. J. Kelly, P. A. Machado, I. Martinez Soler, S. J. Parke, and Y. F. Perez Gonzalez, Phys. Rev. Lett. 123, 081801 (2019), arXiv:1904.02751 [hep-ph] .
  * Kelly  _et al._ (2021) K. J. Kelly, P. A. N. Machado, I. Martinez-Soler, and Y. F. Perez-Gonzalez, (2021), arXiv:2110.00003 [hep-ph] .
  * Kleiss and Pittau (1994) R. Kleiss and R. Pittau, Comput. Phys. Commun. 83, 141 (1994), arXiv:hep-ph/9405257 .
  * Byckling and Kajantie (1969) E. Byckling and K. Kajantie, Nucl. Phys. B 9, 568 (1969).
  * Lepage (1978) G. P. Lepage, J. Comput. Phys. 27, 192 (1978).
  * Lepage (1980) G. P. Lepage, (1980).
  * Strang (1968) G. Strang, SIAM Journal on Numerical Analysis 5, 506 (1968), https://doi.org/10.1137/0705041 .
  * Yoshida (1990) H. Yoshida, Physics Letters A 150, 262 (1990).
  * McLachlan and Quispel (2002) R. I. McLachlan and G. R. W. Quispel, Acta Numerica 11, 341–434 (2002).