# A quantum Monte Carlo based approach to intranuclear cascades

Joshua Isaacson ${}^{\,{\rm a}}$, William I. Jay ${}^{\,{\rm a}}$, Alessandro
Lovato ${}^{\,{\rm b,c}}$, Pedro A. N. Machado ${}^{\,{\rm a}}$, Noemi Rocco
${}^{\,{\rm a,b}}$ ${}^{\,{\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

###### Abstract

We propose a novel approach to intranuclear cascades which takes as input
quantum Monte Carlo nuclear configurations and uses a semi-classical, impact-
parameter based algorithm to model the propagation of protons and neutrons in
the nuclear medium. We successfully compare our simulations to available
proton-carbon scattering data and nuclear-transparency measurements. By
analyzing the dependence of the simulated observables upon the ingredients
entering our intranuclear cascade algorithm, we provide a quantitative
understanding of their impact. Particular emphasis is devoted to the role
played by nuclear correlations, the Pauli exclusion principle, and interaction
probability distributions.

###### pacs:

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

††preprint: FERMILAB-PUB-20-368-T

##  I Introduction

The propagation of nucleons through the nuclear medium is an important aspect
of nuclear reactions, heavy-ion collisions, and astrophysical environments. It
is also crucial in the analysis of electron-nucleus scattering experiments
(see e.g. Refs. Barreau  _et al._ (1983); Makins  _et al._ (1994); Anghinolfi
_et al._ (1995); Abbott  _et al._ (1998); Dutta  _et al._ (2000, 2003); Rohe
_et al._ (2005); Dai  _et al._ (2018, 2019); Murphy  _et al._ (2019)),
neutrino oscillation experiments Abe  _et al._ (2016); Ayres  _et al._ (2007);
Abe  _et al._ (2018); Abi  _et al._ (2020), and dark matter searches Åkesson
_et al._ (2018). For example, scattering of neutrinos on nuclei can produce a
number of outgoing hadrons. The multiplicity of the recoiling hadrons can
indicate statistically if the incoming particle was more likely to be a
neutrino or anti-neutrino Palamara (2016); their momenta are correlated with
the original neutrino energy and direction, which can be used to gather
additional information, for instance, on the leptonic $CP$ phase from the sub-
GeV atmospheric neutrino samples Kelly  _et al._ (2019) or from searches for
dark matter annihilation from the Sun Rott  _et al._ (2017, 2019).

The quantitative understanding of nucleons’ propagation in the nuclear medium
would in principle require a fully quantum-mechanical description of the
hadronic final state. Due to its tremendous difficulty, this problem has been
tackled by introducing different approximations. The seminal papers by Serber
and Metropolis et. al. Serber (1947); Metropolis  _et al._ (1958); Bertini
(1963); Cugnon (1980) laid the foundations for the use of Monte Carlo
techniques in semi-classical intranuclear cascades (INC) that assume classical
propagation between consecutive scatterings. The latter are modeled using
free-space elementary cross sections whose final state is modified to account
for the Pauli principle. The results obtained by these first implementations
of INC agree at least qualitatively with experimental data on nuclear
transparencies and the frequency and angular distribution of emitted fast
protons Metropolis  _et al._ (1958).

Experiments in heavy-ion reactions have spurred theoretical efforts to
describe the dynamical evolution of nucleus-nucleus collisions using transport
methods Bertsch  _et al._ (1984); Stoecker and Greiner (1986); Bauer  _et al._
(1986); Bertsch and Das Gupta (1988); Danielewicz and Bertsch (1991) based on
the Kadanoff-Baym equations Kadanoff and Baym (1962); Botermans and Malfliet
(1990). The impossibility of solving these equations exactly on real-world
computers requires the introduction of approximations, such as the gradient
expansion leading to the Boltzmann-Uehling-Uhlenback equations (BUU). For
instance, the Giessen Boltzmann-Uehling-Uhlenback (GiBUU) model is based on
this truncated set of semi-classical kinetic equations, which describe the
dynamics of the hadronic system explicitly in phase space and in time Cassing
_et al._ (1990); Teis  _et al._ (1997). It can be difficult to estimate
systematic effects coming from the truncation of the transport equations.
Although initially developed to simulate heavy-ion collisions, GiBUU has been
extended to the description of lepton and photon scattering on nuclei Buss
_et al._ (2012); Mosel (2019).

Over the years, several studies have been devoted to improving the accuracy
and extending the predictive power of INC models, with a focus on the analysis
of nuclear spallation processes, where a hadronic probe with energy from a few
tens of MeV to a few hundred MeV strikes a nucleus Cugnon  _et al._ (1997);
Duarte (2007); Iwamoto  _et al._ (2010); Sawada  _et al._ (2012). In contrast
to the standard Glauber approach Glauber  _et al._ (1959), where the so-called
frozen approximation is utilized, INC simulations explicitly account for the
motion of the background particles or scattering centers Golubeva  _et al._
(1998). The validity of the semi-classical propagation assumed in the INC has
yet to be fully assessed. However, genuine quantum-mechanical effects can
either be safely neglected because they are expected to play a minor role—as
for coherent scattering in nuclear transparency calculations—or they can be
effectively parametrized by means of trajectory deflections and nuclear
collective excitations, as in the analysis of nuclear spallation processes
Uozumi  _et al._ (2012).

Important examples of state-of-the-art INC codes used in hadron- and lepton-
nucleus scattering analyses are the Liège INC Boudard  _et al._ (2002), NEUT
Hayato (2002), nuance Casper (2002), PEANUT (used within FLUKA) Battistoni
_et al._ (2013, 2015), NuWro Golan  _et al._ (2012), and GENIE Andreopoulos
_et al._ (2010) programs. Despite some differences in technical aspects and
degree of sophistication, all the above INC models use as input elementary
cross sections and mean-field properties of nuclei, such as single-nucleon
densities. A notable exception is NuWro, which has recently been extended to
incorporate effective nucleon-nucleon correlations Niewczas and Sobczyk
(2019).

In this work, we propose a novel cascade model that employs nuclear
configurations obtained from quantum Monte Carlo (QMC) calculations, which
retains all correlation effects. The probability that the propagating nucleon
scatters with a background particle is modeled using two different weight
functions that depend on the impact parameter of the two nucleons and their
cross section. The nucleon-nucleon cross sections used in the simulation are
taken to be the same as in vacuum; in-medium effects are partially accounted
for by imposing the Pauli exclusion principle and an effective nuclear
binding.

Our semi-classical approach allows for the calculation of exclusive quantities
in nuclear scattering, such as the fully differential phase space and the
number of outgoing nucleons. To test the validity of our model, we compare our
simulations with available data on proton-nucleus scattering cross sections
and nuclear transparencies.

On the technical side, the INC code developed here is publicly available at
https://github.com/jxi24/IntranuclearCascade. The front-end is written in
Python 3, while computationally intensive code is written in C++, using
pybind11 to generate the interface code. The source code is organized to be
portable and modular, making it easy to use, extend, and improve. All
validations are available in the source code and can be easily performed by
the front-end user.

##  II The physics of the impact parameter based intranuclear cascade model

We now discuss in detail the physical content of our INC and how several
ingredients were incorporated in the model. The INC starts by generating a
nuclear configuration, which describes the spatial distribution of protons and
neutrons inside the nucleus. A large number of these configurations are
computed from quantum Monte Carlo (QMC) methods Carlson  _et al._ (2015) that
use as input the highly-realistic Argonne $v_{18}$ (AV18) Wiringa  _et al._
(1995) plus Illinois 7 (IL7) Pieper (2008) Hamiltonian – more details can be
found in Sec. II.1. All particles in a nuclear configuration are initially
labeled as background particles. After randomly picking a nuclear
configuration, a random nucleon is selected, struck, and labeled as a
propagating nucleon. Propagating nucleons are assumed to be point-like, on-
shell particles moving at some velocity given by their four-momentum. The
spatial coordinates of background nucleons are kept fixed until they interact
with a propagating particle. In that case, their momentum is sampled from
either a local or global Fermi gas distribution.

The system is evolved in time steps, which is a parameter of the cascade
model, and should be chosen small enough to simulate the cascade accurately
but large enough to run in a reasonable amount of time. At each time step, the
following procedure is performed for each propagating nucleon. First, it is
checked if there are any background particles within the volume perpendicular
to the line segment defined by the initial to the final position of a
propagating particle; this is schematically shown in Fig. 1.

Figure 1: Schematics of a propagating nucleon (black circle). The distance
travelled by the propagating nucleon is depicted by a black arrow. The crosses
represent background nucleons. An interaction test would be performed for the
nucleon in red.

The impact parameter is calculated for all such nucleons, and this list is
sorted from closest to furthest. Then an accept-reject step is performed on
this list, until either an interaction takes place or the end of the list is
reached. This test determines if an interaction occurred according to a
probability distribution. If an interaction occurred, the phase space is
generated using a fully differential nucleon-nucleon cross section. The Pauli
exclusion principle is approximately taken into account by comparing the
magnitude of the momentum of the final state particles with the Fermi momentum
$k_{F}$ (the user can choose between either the global or the local values).
The interaction takes place only if both momenta are above $k_{F}$. Otherwise,
the interaction does not happen, and the propagating particle keeps its
original four momentum.

If the interaction took place, the outgoing particles are both treated as
propagating particles, and a formation zone is set for each of them Landau and
Pomeranchuk (1953); Stodolsky (1975) (see also Ref. Golan  _et al._ (2012)).
Note that, as a distinctive feature of our INC model, the interaction is
finite ranged and the outgoing nucleons do not have to be at the same position
in space. Finally, a propagating particle that reaches a sufficiently large
radius may exit the nucleus if its kinetic energy overcomes the effective
nuclear binding. If there is insufficient energy to escape, the particle is
relabeled as a background particle. Otherwise, this effective binding energy
is subtracted from the particle’s energy, and then the particle is labeled as
a final particle and stops propagating. The INC stops when there are no more
propagating particles in the nucleus.

The structure of the algorithm is summarized in Fig. 2. In what follows, we
provide the details and expressions that are used in our impact parameter
based INC.

Start Generate nuclear configuration Kick random nucleon Evolve propagating
nucleon(s) for time $\delta t$ Possible interaction? Outside nucleus? Remove
from propagating particle list Generate outgoing momenta Pauli Blocked?
Restore original outgoing momenta Update outgoing momenta of particles Add
particle to list and update formation zone NoYesYesNoYesNo Figure 2:  The
proposed algorithm for the INC model.

###  II.1 Nuclear configuration

Nuclei are complicated many-body systems of fermions, whose structure and
dynamics emerge from individual interactions among the constituent protons and
neutrons. To a remarkably large extent, the latter can be modeled by the non-
relativistic Hamiltonian

$H=\sum_{i}\frac{{\bf
p}_{i}^{2}}{2m_{N}}+\sum_{i<j}v_{ij}+\sum_{i<j<k}V_{ijk}\,,$ (1)

where ${\bf p}_{i}$ denotes the momentum of the $i$-th nucleon with mass
$m_{N}$, while $v_{ij}$ and $V_{ijk}$ are the nucleon-nucleon ($N\\!N$) and
three-nucleon ($3N$) potentials, respectively. We employ the highly-realistic
Argonne $v_{18}$ (AV18) $N\\!N$ interaction Wiringa  _et al._ (1995) that
includes spin, isospin, tensor, spin-orbit and quadratic momentum-dependent
terms as well as isospin-symmetry-breaking corrections and reproduces the
Nijmegen nucleon-nucleon database with a $\chi^{2}$/datum $\simeq 1$. For the
$3N$ interaction, we use the Illinois-7 (IL7) potential Pieper (2008), which
consists of the dominant standard Fujita-Miyazawa two-pion exchange and
smaller multi-pion-exchange components resulting from the excitation of
intermediate $\Delta$ resonances. The IL7 potential also contains
phenomenological isospin-dependent central terms. The parameters
characterizing this three-body potential have been determined by fitting the
low-lying spectra of nuclei in the mass range $A$=3–10. The overall AV18+IL7
Hamiltonian then leads to predictions of $\approx 100$ ground- and excited-
state energies up to $A$=12, including the 12C ground- and Hoyle-state
energies, in good agreement with the corresponding empirical values Carlson
_et al._ (2015).

Over the last decades significant progress has been made in the development of
ab initio nuclear methods, which solve the many-body Schrödinger equation
associated with the Hamiltonian of Eq. (1) with controlled approximations
Bogner  _et al._ (2010); Barrett  _et al._ (2013); Hagen  _et al._ (2014);
Barbieri and Carbone (2017). For light nuclei, quantum Monte Carlo (QMC) and,
in particular, Green’s function Monte Carlo (GFMC) methods have been exploited
to carry out calculations of nuclear properties, based on realistic $N\\!N$
and $3N$ potentials, and consistent one- and two-body meson-exchange currents
Carlson  _et al._ (2015). GFMC begins with the construction of a trial wave
function $\Psi_{T}$ that is a symmetrized product of two- and three-body
correlation operators acting on an antisymmetric $A$-body single-particle wave
function that has the proper quantum numbers for the state of interest. The
variational parameters in $\Psi_{T}$ are found by minimizing the energy
expectation value

$E_{0}\leq
E_{T}=\frac{\langle\Psi_{T}|H|\Psi_{T}\rangle}{\langle\Psi_{T}|\Psi_{T}\rangle}\,,$
(2)

where $E_{0}$ is the true ground-state energy of the system. The calculation
of $E_{T}$ requires the numerical solution of a multidimensional integral that
is carried out employing standard Metropolis Monte Carlo sampling in
configuration space.

GFMC then projects out the lowest eigenstate $\Psi_{0}$ of the given quantum
numbers starting from $\Psi_{T}$ by performing a propagation in imaginary time
$\tau$

$|\Psi_{0}\rangle=\lim_{\tau\rightarrow\infty}{\rm
exp}[-(H-E_{0})\tau]|\Psi_{T}\rangle.$ (3)

The propagation $|\Psi(\tau)\rangle={\rm exp}[-(H-E_{0})\tau]|\Psi_{T}\rangle$
is carried out as a series of many small imaginary-time steps $\Delta\tau$.
Expectation values of operators are evaluated as mixed matrix elements
$O(\tau)=\langle\Psi_{T}|O|\Psi(\tau)\rangle$, and the behavior as a function
of $\tau$ analyzed to obtain converged results. Because $H$ and ${\rm
exp}[-(H-E_{0})\tau]$ commute, the mixed estimate is the exact expectation of
$\langle\Psi(\tau/2)|O|\Psi(\tau/2)\rangle$ but linear extrapolations are used
to evaluate other quantities.

In addition to binding energies the GFMC provides detailed information on the
distribution of nucleons in a nucleus in both coordinate and momentum space,
which are interesting in multiple experimental settings. For example, the
mixed-estimate of the single-nucleon density is calculated as

$\displaystyle\rho_{N}(r)$ $\displaystyle=\frac{1}{4\pi
r^{2}}\big{\langle}\Psi_{T}\big{|}\sum_{i}\mathcal{\delta}(r-|\mathbf{r}_{i}|)P_{N}\big{|}\Psi(\tau)\big{\rangle}\,,$
(4)

where $N=p,n$; $P_{N_{i}}=\frac{1\pm\tau_{z_{i}}}{2}$ is the neutron or proton
projector operator; and, $\rho_{N}$ integrates to the number of protons or
neutrons. The two-body density distribution, yielding the probability of
finding two nucleons with separation $r$, is defined as

$\displaystyle\rho_{NN}(r)$ $\displaystyle=\frac{1}{4\pi
r^{2}}\big{\langle}\Psi_{T}\big{|}\sum_{i<j}\mathcal{\delta}(r-|\mathbf{r}_{ij}|)P_{N_{i}}P_{N_{j}}\big{|}\Psi(\tau)\big{\rangle}\,.$
(5)

The positions of the constituents protons and neutrons utilized in the nuclear
cascade algorithm are sampled from $36000$ GFMC configurations. We employ the
so-called constrained-path approximation Wiringa  _et al._ (2000) to make sure
that their Monte Carlo weights remain positive, thereby facilitating their
usage in the cascade algorithm. As a consequence, the single-proton
distribution displayed by the blue solid circles of Fig. 3 is slightly
different from the results reported in Ref. Lovato  _et al._ (2013), which
have been obtained performing fully unconstrained imaginary-time propagations.
Since we neglect the charge-symmetry breaking terms in the Hamiltonian, and
since 12C is isospin symmetric, the single-neutron distribution is identical
to that of the proton.

For benchmark purposes, we also sample $36000$ mean-field (MF) configurations
from the single-proton distribution. The corresponding single-proton densities
coincide by construction with the GFMC one, as shown in Fig. 3. However, the
differences between GFMC and MF configurations become apparent when comparing
the corresponding two-body density distributions represented in Fig. 4. The
short-range repulsive core of the $NN$ interaction prevents two nucleons from
being close to each other. As a consequence, the $pp$ and $np$ GFMC density
distributions are small at short separation distances. Furthermore, the
difference between the GFMC $pp$ and $np$ density distributions around $r=1$
fm can be attributed to the strong tensor correlations induced by the one-
pion-exchange part of the $NN$ interaction, which is further enhanced by the
two-pion-exchange part of the $3N$ potential. Note that the short-range
behavior of $\rho_{NN}$, which is largely nucleus independent, does depend
strongly on the $NN$ interaction model Cruz-Torres  _et al._ (2019). On the
other hand, the MF ones do not exhibit this rich behavior as the correlations
among nucleons are entirely disregarded.

Figure 3: Nucleon density in carbon from Green’s function Monte Carlo (red)
and mean field (blue) configurations.

Figure 4: Proton-proton (top panel) and proton-neutron (bottom panel)
correlation functions in carbon from Green’s function Monte Carlo (red) and
mean field (blue) configurations.

###  II.2 Nucleon momentum distribution

As mentioned above, when a nucleon is struck, its momentum is obtained
assuming either a local or global Fermi gas distribution. In the case of the
local Fermi gas, the magnitude of the three-momentum is randomly sampled in
the interval $[0,k^{N}_{F}(r)]$ where $k^{N}_{F}(r)$ is the Fermi Momentum
defined in terms of the single nucleon density
$k^{N}_{F}(r)=(\rho_{N}(r)3\pi^{3})^{1/3}$ and $N=p,n$. In the case of the
global Fermi gas, the momentum is determined in the same way, but $k^{N}_{F}$
is position independent. The local Fermi gas model is known to provide a more
realistic nucleon momentum distribution for finite nuclei than the global
Fermi gas. For this reason, although both models are implemented in our code,
we only present results based on the local Fermi gas predictions. In the
future, we plan to include more accurate nucleon momentum distribution, based
on state-of-the-art many-body calculations that properly account for nuclear
correlations.

###  II.3 Nucleon-nucleon interaction algorithm

To check if an interaction between nucleons occurs, an accept-reject test is
performed on the closest nucleon according to a probability distribution
$P(b)$ (see e.g. Ref. Sjostrand and van Zijl (1987) for similar
considerations) where $b$ is the impact parameter. We impose two conditions on
this probability,

$P(0)=1\quad\textrm{and}\quad\int_{0}^{2\pi}\int_{0}^{\infty}d\varphi\,bdbP(b)=\sigma,$
(6)

where the cross section $\sigma$ depends on the incoming particle content and
the center-of-mass energy, which is sampled from the nuclear configuration.
The second condition ensures that the mean free path of a nucleon traveling in
a medium of uniform density is $\lambda_{\rm mfp}=1/\sigma\bar{\rho}$, where
$\bar{\rho}$ is the number density.

Two implementations of $P(b)$ have been studied here. The first we dub the
_cylinder_ interaction probability,

$P_{\rm cyl}(b)=\Theta(\sigma/\pi-b^{2}),$ (7)

where $\Theta(x)=1$ if $x\geq 0$, else $\Theta(x)=0$. This probability mimics
a more classical, billiard ball like system, where each billiard ball has a
radius $\approx\sqrt{\sigma/\pi}$. The second implementation is the _Gaussian_
interaction probability

$P_{\rm Gau}(b)\equiv\exp\left(-\frac{\pi b^{2}}{\sigma}\right),$ (8)

which is inspired by the work of Ref. Sjostrand and van Zijl (1987). Both
$P_{\rm cyl}$ and $P_{\rm Gau}$ satisfy the conditions in Eq. (6). We use the
nucleon-nucleon cross sections from the SAID database SAID database  obtained
using GEANT4 Agostinelli  _et al._ (2003), or from the NASA parametrization
Norbury (2008).

###  II.4 Phase space, Pauli blocking and after-interaction

If an interaction occurred, the phase space of the outgoing particles is
generated using fully differential nucleon-nucleon cross sections. Note that,
at the moment, we only include protons and neutrons in our INC model. Pauli
blocking enforces Fermi-Dirac statistics for the nucleons and amounts to
testing whether their final-state momenta are above the Fermi momentum. Two
different models of the Pauli exclusion principle have been approximately
implemented. The _global_ and _local_ Pauli blocking routines essentially
forbid a scattering if the momentum of any of the final state particles is
below the average Fermi momentum (for the global Fermi gas model) or the local
Fermi momentum (for the local Fermi gas model), respectively. We emphasize
again that, although we have implemented the global Fermi gas model, we do not
report any results using it.

If the interaction took place, the outgoing particles are both treated as
propagating particles, and a formation zone is set for each of them Landau and
Pomeranchuk (1953); Stodolsky (1975) (see also Ref. Golan  _et al._ (2012)).
The formation zone is a length or period of time in which a particle does not
interact with any nucleons. This models the coherence and interference of
interactions in quantum mechanics. The formation zone is given by

$\delta t=\frac{E^{\prime}}{m_{N}^{2}-p\cdot p^{\prime}},$ (9)

where $m_{N}=938$ MeV is the nucleon mass, $E^{\prime}$ is the energy of the
outgoing nucleon, and $p$ and $p^{\prime}$ are the four momenta of the
incoming and outgoing nucleon respectively.

###  II.5 Exiting the nucleus

When the radial position of a propagating particle is larger than a certain
distance, much larger than the nuclear radius, a test is performed to check if
the particle has enough energy to escape the nuclear potential. If its kinetic
energy is larger than the the nuclear potential barrier, then the particle
escapes the nucleus and is labeled as a final state particle. The momentum of
the particle is modified to include the effective nuclear binding. Final state
particles do not propagate and are essentially the input that should be given
to a detector simulation.

##  III Results and Model Validation

In this Section, we present the validation tests we performed and the results
we obtained within our INC model. For comparison purposes, we also implement a
version of the nucleon-nucleon interaction algorithm that we dub _mean free
path_ (MFP). This approach is routinely used in event generators. Within this
algorithm, the system is not evolved in time steps but rather in constant
position steps that we indicate as $\lambda_{\rm max}$. At each step of the
loop the mean free path of the propagating particle reads

$\tilde{\lambda}=1/\Big{(}\rho_{p}(r)\sigma_{Np}+\rho_{n}(r)\sigma_{Nn}\Big{)}$
(10)

where $N=n,p$ refers to the isospin of the propagating particle and $r$ its
distance from the centre of the nucleus. Note that in the limit of constant
density, $\rho_{p(n)}(r)=\bar{\rho}$, we would recover
$\tilde{\lambda}=\lambda_{\rm mfp}$. (see also Sec. III.1 below). The
probability that the struck particle traveled a distance $\lambda$ without
interacting can be written as $P(\lambda)=e^{-\lambda/\tilde{\lambda}}$. We
can then sample $\lambda=-\tilde{\lambda}\cdot\ln({\rm rand}[0,1])$ and say
that an interaction took place if $\lambda<\lambda_{\rm max}$. Note that
$\lambda_{\rm max}$ has to be chosen small enough in order to satisfy the
assumption of a constant density.

Figure 5: Distribution of distance traveled (blue histogram) by 500 MeV test
nucleons traveling through a $0.16$ nucleons/fm3 background density of
nucleons at rest, fixed interaction cross-section of 50 mb, and using the
Gaussian interaction probability. The black and orange lines are the expected
distributions for the fitted (1.23 fm) and expected (1.25 fm) values of the
mean free path, respectively.

###  III.1 Mean free path

The simplest validation consists of calculating the mean free path of a
nucleon travelling through a medium of uniformly distributed nucleons, with
fixed interaction cross sections. Our goal is to obtain

$\lambda_{\rm mfp}=\frac{1}{\sigma\rho},$ (11)

where $\sigma$ is the nucleon-nucleon cross section, and $\rho$ is the number
density of nucleons. We simulate a background of static, randomly and
uniformly distributed protons and neutrons at rest. For an event, we insert
into the medium a test particle with a fixed energy and propagate it through
the medium. When the first interaction happens, we record the length the test
particle has traveled. Figure 5 shows the distribution of lengths traveled for
a 500 MeV test proton in a $0.16$ nucleons/fm3 background density of nucleons,
assuming a fixed interaction cross section of 50 mb, and using the Gaussian
interaction probability. The expected distributions, according to Eq. (11),
are also shown for the fitted (black line) and the expected (orange line)
values of the mean free path. The nuclear density and interaction cross
section were chosen arbitrarily. Adjusting these values does not change the
agreement between the expected mean free path and the calculated mean free
path. While not shown here, using the cylinder interaction probability does
not change the results either. As we can see, our code correctly reproduces
the expected behavior for the mean free path, allowing us to proceed to more
complex tests of our INC.

###  III.2 Proton-carbon Scattering Data

Reproducing the proton-nucleus cross section measurements is an important test
of the accuracy of the INC model. Proton-nucleus scattering probes the
nucleon-nucleon cross section which is typically divided into two pieces, the
reaction and the elastic cross sections,

$\sigma_{\rm tot}=\sigma_{R}+\sigma_{\rm el}.$ (12)

In the elastic part, no energy is transferred into nuclear excitation and the
nucleus remains unbroken, that is $n+A\rightarrow n+A$. The reaction cross
section includes transition to nuclear excited states, $n+A\rightarrow
n+A^{*}$, as well as inelastic reactions $n+A\rightarrow X$.

Several experiments have been carried out to determine the total reaction
cross section, see for example Refs. Dicello and Igo (1970); Bobchenko  _et
al._ (1979); Bauhoff (1986); Kox  _et al._ (1987); Sihver  _et al._ (1993);
Carlson (1996). The latter is typically obtained by measuring the total cross
section from the change in intensity of a calibrated proton beam traversing a
carbon target and then subtracting the calculated elastic cross section.

Figure 6: Proton-carbon scattering total cross section as a function of the
incoming proton kinetic energy. In the upper panel the entire energy range for
which experimental data are available is shown. In the lower panel the low
energy region is magnified. The red and blue curves correspond to the cylinder
algorithm where the mean field (MF) and quantum Monte Carlo (QMC)
configurations have been used, respectively. The green and orange curves are
the same but for the Gaussian interaction probability. The results displayed
in purple refers to the mean free path (MFP) calculations. The solid and
dashed curves corresponds to the use of the GEANT4 Agostinelli  _et al._
(2003) and NASA Norbury (2008) parametrization of the cross section in the
interaction probability, respectively. The data points are from Ref. Carlson
(1996)

We compute $\sigma_{R}$ neglecting Coulomb interactions, as they are expected
to contribute mostly to $\sigma_{\rm el}$. We obtain the proton-carbon
scattering cross section by the following simulation (with a different setup
from the proposed algorithm of Fig. 2). We define a beam of protons with
energy $E$, uniformly distributed over an area $A$ (orthogonal to the proton
momenta). Note that $A\gg\pi R^{2}$, where $R$ is the radius of the carbon
nucleus. The carbon nucleus is situated in the center of the beam. We
propagate each proton in time and check for scattering at each step. The Monte
Carlo reaction cross section is then defined as the area of the beam times the
fraction of scattered events, namely,

$\sigma_{\rm MC}=A\frac{N_{\rm scat}}{N_{\rm tot}}.$ (13)

This is not exactly the experimentally measured reaction cross section.
Angular and/or momentum acceptances for the attenuated beam are finite, and we
do not include these effects in our calculation. Nevertheless, we do not
expect such effects to change our results significantly, and thus $\sigma_{\rm
MC}$ should be a good approximation of the reaction cross section. Moreover,
imposing Pauli blocking on both outgoing nucleons will effectively suppress
the contribution of elastic transitions.

The two panels of Fig. 6 display the proton-carbon scattering cross sections
as a function of the proton kinetic energy. In the upper panel our Monte Carlo
simulations are compared with experimental data in the entire energy region in
which data are available Carlson (1996), while the lower panel focuses on
proton kinetic energies below 200 MeV. The curves correspond to different
implementations of the INC. These implementations are composed of three
ingredients, namely,

  1. 1.

Nuclear configuration: quantum Monte Carlo (QMC) or mean field (MF);

  2. 2.

Interaction model: cylinder (cyl), Gaussian (Gauss), or mean free path (MFP);

  3. 3.

Nucleon-nucleon cross section: elastic (El) or total (Tot).

Note that the interaction model MFP does not use any configuration, but rather
the local density.

The solid lines have been obtained using the nucleon-nucleon cross sections
from the SAID database SAID database , obtained using GEANT4 Agostinelli  _et
al._ (2003), in which only the elastic contribution is retained. The dashed
lines used the NASA parameterization Norbury (2008), which includes
inelasticities. The inelastic contribution leads to an enhancement of the
cross section which is necessary to reproduce the data correctly for
intermediate to large proton kinetic energies. The NASA parameterization
allows us to approximate the effects of including pions in the cascade.
Results obtained with these two nucleon-nucleon cross sections are indicated
in Fig. 6 by “El” for elastic cross section and “Tot” for the total cross
section. As expected, the $p$-carbon cross section obtained using only the
elastic nucleon-nucleon cross section is consistently lower than the one
obtained using the total nucleon-nucleon cross section. The effect is large
for $T_{p}>50$ MeV, where pion production becomes relevant. In a future work,
we will explicitly include pion degrees of freedom in our INC as these are
known to play a crucial role in this region. Additionally, the present results
neglect modifications to the nucleon-nucleon cross section arising from the
nuclear medium itself. Such modifications have been considered in the
literature Pandharipande and Pieper (1992). We also leave these considerations
for a future work.

The dependence of the results on the functional form of the interaction
probability and on the nuclear model adopted to generate nuclear
configurations has also been investigated. The blue and red lines are obtained
from the cylinder probability of Eq. (7) and sampling the initial nucleon
configurations from the quantum Monte Carlo (QMC) and mean field (MF)
distributions, respectively. To obtain the orange and green curves, QMC and MF
configurations were used together with the Gaussian probability of Eq. (8). We
observe that in both cases, the QMC and MF curves are almost superimposed,
indicating that this observable does not depend strongly on correlation
effects among the nucleons. A more detailed discussion of their impact on
particle propagation distances will be discussed in Sec. III.4.

The purple lines refer to mean free path (MFP) calculations which present
conceptual differences with respect to the cylinder and Gaussian cases and
yield large discrepancies in the results. The predictions for the $p$-carbon
cross section obtained from the MFP approach underestimate the experimental
data and are significantly lower with respect to the other curves obtained
with more refined techniques.

In the lower panel of Fig. 6 we focus on the comparison between the different
MC simulations and the data in the low energy region, i.e., proton kinetic
energies below $200$ MeV. The results obtained using the Gaussian and cylinder
algorithm are in fair agreement with the experimental data: the curves display
the correct behaviour although the position of the peak is not exactly
reproduced. We note that Pauli blocking plays a fundamental role for these
kinematics, significantly quenching the results. More sophisticated techniques
aimed at consistently orthogonalizing the nuclear wave function in the final
state are discussed in Ref. González-Jiménez  _et al._ (2019). Their
implementation in our INC will be the subject of a future work.

The functional forms adopted to determine the probability distribution
(Gaussian vs. cylinder) yield predictions which differ by $\sim 15\%$ in the
low energy region (compare, for instance, the green and red lines in Fig. 6).
For intermediate proton kinetic energies $T_{p}$ the curves converge to the
same result and reproduce the data. At low $T_{p}$, the predicted cross
sections for the MFP method are $\sim 50\%$ smaller than those of the Gaussian
and cylinder methods. As these differences are also present in nuclear
transparency, we will discuss their origins in the next section.

###  III.3 Nuclear transparency

The nuclear transparency yields the average probability that a struck nucleon
(either a proton or a neutron) leaves the nucleus without interacting with the
spectator particles. Measurements of the nuclear transparency to high energy
protons in quasielastic $e+A\to e+p+A^{\prime}$ scattering have been carried
out by a number of experiments O’Neill  _et al._ (1995); Abbott  _et al._
(1998); Rohe  _et al._ (2005); Garrow  _et al._ (2002); Dutta  _et al._
(2003); Garino  _et al._ (1992). These measurements are performed with fixed
kinematics (i.e., for a fixed incoming beam energy and scattering angle) and
measure the outgoing electron and proton. The proton produced in the primary
vertex can be absorbed or deflected while exiting the nuclear medium because
of final state interactions with the remnant system, leading to a reduction of
the measured $e+A\to e+p+A^{\prime}$ cross section. This reduction defines the
nuclear transparency, which is given by the ratio of the observed events to
events predicted in the Plane Wave Impulse Approximation (PWIA)

$T=\frac{N_{\rm exp}}{N_{\rm PWIA}}.$ (14)

To obtain $N_{\rm PWIA}$ the initial proton is treated as a bound state with
energy and momentum distribution described by a spectral function; while the
final proton is a free particle state propagating as a plane wave.

Figure 7: Carbon transparency as a function of the proton kinetic energy. The
different curves indicate different approaches used as described in Fig. 6.
The experimental data are taken from Refs. O’Neill  _et al._ (1995); Abbott
_et al._ (1998); Rohe  _et al._ (2005); Garrow  _et al._ (2002); Dutta  _et
al._ (2003); Garino  _et al._ (1992)

In our Monte Carlo simulation, the nuclear transparency has been obtained via
the following procedure. We randomly sample a nucleon, its position being
generated according to either the QMC or MF distributions. We give the nucleon
a kick by assigning it a given kinetic energy $T_{p}$ and three-momentum ${\bf
p}$, and propagate it through the nuclear medium. The Monte Carlo transparency
is then defined as

$T_{\rm MC}=1-\frac{N_{\rm hits}}{N_{\rm tot}}$ (15)

where Pauli blocking has been implemented in the determination of $N_{\rm
hits}$. Note that for a given initial and final energy and scattering angle of
the electron, one can unambiguously define the momentum ${\bf q}$ transferred
to the target nucleus. The direction and the momentum of the nucleon in the
final state has to be determined applying energy- and momentum-conservation
relations and accounting for the Fermi motion of the struck nucleon in the
initial state. It follows that defining the kinematics of the hadronic final
state after the hard scattering depends on the nuclear model of choice.
However, in the analysis of different experiments, the data are given as a
function of the average nucleon momentum (and kinetic energy) given by ${\bf
p}={\bf q}$ ($T_{p}=\sqrt{|{\bf q}|^{2}+m_{N}^{2}}-m_{N}$).

In Fig. 7 we compare the nuclear transparency data from Refs. O’Neill  _et
al._ (1995); Abbott  _et al._ (1998) to our predictions. The different lines
are the same as for Fig. 6. We find an overall satisfactory agreement between
the Gaussian and cylinder curves with the experimental data once inelastic
effects are taken into consideration; this corresponds to the results using
the NASA parametrization for the nucleon-nucleon cross sections. For moderate
to large values of the proton kinetic energy, pions play an important role in
quenching the transparency. Moreover, the Gaussian and cylinder curves exhibit
correct behavior consistent with the data also for small $T_{p}$ where the
simplified MFP model described above fails. As in Fig. 6, we observe very
small differences between the QMC and MF calculations. For low and
intermediate kinetic energies, the transparency obtained from the MFP approach
is much smaller than the corresponding results for the cylinder and Gaussian
curves.

Figure 8: Left panel: a schematic picture of an external proton scattering off
the nucleus. The distance from the proton to the center of the nucleus is
$r_{1}$, and the propagation step is $d\ell$. The radius of the cylinder is
given by $\sqrt{\sigma/\pi}$ where $\sigma$ is the interaction cross section
between the proton and a background particle; $d\ell$ is also the height of
the cylinder. Right panel: same as for the left one, but for a nucleon kicked
inside the nucleus. This follows what is done in the nuclear transparency
event simulations.

Finally, we discuss the origin of the discrepancies between the MFP and the
cylinder algorithm with MF configurations for the $p$-carbon cross section and
carbon transparency. Both approaches rely on the single-nucleon density
distribution to sample the initial nucleon positions (nuclear correlations are
neglected) but use different definitions of the interaction probability. The
left panel of Fig. 8 schematically shows one contribution to the $p$-carbon
cross section in which the proton is at a distance $r_{1}$ larger than the
nuclear radius. In the cylinder algorithm, the interaction probability is
equal to one if a particle is present in the volume defined by:
$V=d\ell\cdot\sigma$. Both $\sigma_{pp}$ and $\sigma_{np}$ have a maximum for
low proton momentum values. Hence, for low momenta, the probability of
interaction could be non-vanishing even when the projectile proton is far from
the center of the nucleus.
On the other hand, within the MFP approach, if the probe is outside the
nucleus then the approximation of a constant density $\rho(r_{1})=0$ within
the volume $V=d\ell\cdot\sigma$ yields a vanishing interaction probability.
This different behaviour leads to a lower $p$-carbon cross section using the
MFP approach, as observed in Fig. 6. When computing the nuclear transparency
we kick a nucleon which is located inside the nucleus as displayed in the
right panel of Fig. 8. In this case, assuming a constant density is more
likely to overestimate the interaction probability, especially for low momenta
where the cross section is larger. This observation is consistent with Fig. 7
where the MFP curves predict a larger number of interactions, and therefore a
lower nuclear transparency, for small $T_{p}$.

###  III.4 Correlation effects

Figure 9: The four panels corresponds to histograms of the distance traveled
by a struck particle before the first interaction takes place for different
values of the interaction cross section. The results in blue and red
correspond to MF and QMC initial nucleon configurations, respectively. For
each of the panels we also report the fixed cross-section used, the total
number of events generated, and the number of hits for each configuration.

The role played by nuclear correlations in final state interactions of the
recoiling nucleon has been investigated in Refs. Pandharipande and Pieper
(1992); Radici  _et al._ (1994); Benhar  _et al._ (1993); Nikolaev  _et al._
(1993); Benhar (2013). As discussed in Ref. Benhar  _et al._ (1995) the hit
nucleon is surrounded by a short-distance correlation hole produced by both
the Pauli principle and the repulsive nature of realistic nuclear
interactions. Because of this correlation hole, the stuck nucleon is expected
to freely propagate for $\sim 1$ fm before interacting with any of the
background particles. To test the validity of these observations in our INC
model, in Fig. 9 we report the histograms of the distance traveled by a struck
nucleon before its first interaction occurs—we stop the simulation
afterwards—with each panel corresponding to a different value of the
interaction cross section. In order to gauge the effect of nuclear
correlations, the initial positions of the nucleons are sampled from either MF
(blue) or QMC (red) configurations. A random nucleon inside the nucleus is
recoiled and assigned a momentum of 200 MeV. Pauli Blocking has been neglected
here to isolate the dependence of the results on the spatial distribution of
the nucleons. We employ the cylinder algorithm and use a fixed cross
section—which determines the cylinder base area—varying between 0.5 and 100
mb.

For $\sigma=0.5$ and 10 mb, the volume spanned by the propagating particle is
very small. The first and second panels of Fig. 9 clearly show the MF
distribution peaking toward smaller distances than the QMC distribution. This
difference primarily originates from the short-range repulsion of the AV18
potential that reduces the probability of finding two nucleons close to each
other and allows the struck particle to propagate longer before interacting.
This effect is more pronounced for cross sections below about $10~{}{\rm
mb}=1~{}{\rm fm}^{2}$ since correlations affect nucleon configuration for
inter-particle distances within $1\sim 2$ fm, as can be seen in Fig. 4. On the
other hand, larger cross sections yield larger cylinders. In this case, the
propagating particle becomes less sensitive to the local distribution of
nucleons and more sensitive to the integrated density in a larger volume,
reducing the effect of correlations. For these larger cross sections, the MF
and QMC event distributions follow the same trend, as can be seen in the lower
panels of Fig. 9, corresponding to $\sigma=50$ and $100$ mb.

In each panel we also report the number of hits and the total number of events
simulated. Obviously, the number of interactions increases for the larger
cross sections. It is more interesting to notice that the total number of hits
does not present a strong dependence on configurations with (QMC) and without
(MF) correlations. This explains why sampling from these two different set of
configurations give comparable results for the $p$-carbon cross section and
carbon transparency, displayed in Figs. 6 and 7. These observables are indeed
more sensitive to the occurrence of any interaction rather than to its
location inside the nucleus. This is especially true for the $p$-carbon cross
section, since the attenuation of a proton beam is sensitive to the total
nuclear density rather than its differential profile (compare, e.g., the solid
blue and red lines in Fig. 6).

An alternative way to include correlation effects in INC models, based on two-
body density distributions, has been recently proposed by the NuWro
collaboration and leads to a larger enhancement of the nuclear transparency
than we find Niewczas and Sobczyk (2019). This difference can partly be
ascribed to correlation terms between two spectator nucleons Nikolaev  _et
al._ (1993) that are automatically incorporated in our QMC configurations, but
cannot be modeled by the two-body correlation functions alone. In addition, as
discussed in Sec. III.4, when the nucleon-nucleon cross section is larger than
about $20$ mb, the base of the propagating cylinder in our model covers the
correlation hole and reduces the importance of correlations. This effect is
absent when only local properties of two-body distribution functions are
considered.

Although nuclear correlations in final state interactions do not play a
prominent role in the observables that we have considered, correlations may be
important in other experimentally relevant quantities or in high-precision
studies. These observables include proton multiplicity and the distribution of
outgoing protons’ direction and energy. We plan to continue exploring these
questions in future work.

##  IV Conclusions

We have presented a novel INC model that takes as input realistic QMC nuclear
configurations to sample the distributions of protons and neutrons inside the
nucleus. Either a cylinder or a Gaussian distribution is used to define the
interaction probability of the propagating nucleon as a function of the impact
parameter and the nucleon-nucleon cross sections. We have considered the
elastic and the total cross sections, corresponding to the SAID database SAID
database  and to the NASA parametrization Norbury (2008), respectively.

To validate the INC we have compared our simulations with experimental data
for $p$-carbon scattering and carbon transparency. For both observables, we
find minimal difference between the cylinder and Gaussian distributions, and
both correctly reproduce the trend of the data. We gauge the role of nuclear
correlations in the propagation of the recoiling nucleons by utilizing nuclear
configurations computed within quantum Monte Carlo and comparing with mean
field distributions. The distance traveled by a struck particle before the
first interaction occurs, for sufficiently small value of the nucleon-nucleon
cross section, is significantly longer when nuclear correlations are accounted
for. This is due to the short-range repulsion characterizing the AV18
potential that gives rise to a correlation hole surrounding the propagating
nucleon, thereby reducing the probability of finding two nucleons close to
each other. However, the total number of hits is only mildly affected by the
presence of nuclear correlations. As a consequence, they mildly affect our
calculations of $p$-carbon scattering and carbon transparency. The results
presented for the latter quantity by the NuWro collaboration indicate a larger
importance of nuclear correlations, whose implementation is based on the two-
body density distributions Niewczas and Sobczyk (2019).

Nucleon-nucleon cross sections are an important ingredient in INC models.
Since the scattering between the two nucleons takes place in the nuclear
medium, the free nucleon-nucleon cross sections have to be corrected. In our
simulation, we do so by approximately implementing the Pauli exclusion
principle: in the aftermath of each collision, the magnitude of final
nucleons’ momenta have to be larger than the Fermi momentum. This constraint
reduces the effective cross section and considerably improves the agreement
between our simulations and experimental data. As a next step, along the line
of Ref. Pandharipande and Pieper (1992), we plan to replace the bare nucleon
mass with an effective one that depends on the momentum and the nuclear
density, leading to a more complete and consistent treatment of in-medium
effects.

At energies larger than the pion-production threshold, inelastic contributions
to the nucleon-nucleon cross section become relevant. Using the NASA
parametrization, which includes these inelasticities, noticeably improves the
agreement with $p$-carbon scattering and carbon transparency experimental data
in the higher-energy region. Our INC model has yet to include pion propagation
in the nuclear medium. Existing event generators have to rely on a number of
semi-phenomenological approaches constrained by experimental data Salcedo  _et
al._ (1988); Buss  _et al._ (2012); Ashery  _et al._ (1981, 1984); Jones  _et
al._ (1993), the validity of which strongly depends upon the energy of the
propagating pion. A consistent implementation of pion propagation and
absorption is an extremely challenging problem and we leave this aspect for
future works.

The transparency measurements used here depend upon theoretical calculations
based on the plane wave impulse approximation. To compare directly to the
experimental data, we plan on implementing the primary interaction vertex
using the spectral function formalism Benhar  _et al._ (1994); Rocco  _et al._
(2016); Benhar  _et al._ (2008). Additionally to further validate our model,
we will carry out extensive comparison against available semi-exclusive
electron scattering data Cruz-Torres  _et al._ (2020); Duer  _et al._ (2019).

##  V Acknowledgments

We thank John Arrington, Omar Benhar, Donald F. Geesaman, Stefan Höche, T.S.
Harry Lee, Kajetan Niewczas, Stefan Prestel, Jan Sobczyk, and Mike Wagman for
useful discussions . PM and AL thank the Argonne Physics Division and the
Fermilab theory group, respectively, for providing consistent office space
necessary for the completion of this work. 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. and N.R.) and by the NUCLEI SciDAC program.

## References

  * Barreau  _et al._ (1983) P. Barreau _et al._ , Nucl. Phys. A 402, 515 (1983).
  * Makins  _et al._ (1994) N. Makins _et al._ , Phys. Rev. Lett. 72, 1986 (1994).
  * Anghinolfi  _et al._ (1995) M. Anghinolfi _et al._ , J. Phys. G 21, L9 (1995).
  * 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).
  * 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 .
  * Dai  _et al._ (2018) H. Dai _et al._ (Jefferson Lab Hall A), Phys. Rev. C 98, 014617 (2018), arXiv:1803.01910 [nucl-ex] .
  * Dai  _et al._ (2019) H. Dai _et al._ , Phys. Rev. C 99, 054608 (2019), arXiv:1810.10575 [nucl-ex] .
  * Murphy  _et al._ (2019) M. Murphy _et al._ , Phys. Rev. C 100, 054606 (2019), arXiv:1908.01802 [hep-ex] .
  * Abe  _et al._ (2016) K. Abe _et al._ (T2K), Phys. Rev. D93, 112012 (2016), arXiv:1602.03652 [hep-ex] .
  * Ayres  _et al._ (2007) D. S. Ayres _et al._ (NOvA),  (2007), 10.2172/935497.
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE), (2020), arXiv:2002.02967 [physics.ins-det] .
  * Åkesson  _et al._ (2018) T. Åkesson _et al._ (LDMX), (2018), arXiv:1808.05219 [hep-ex] .
  * Palamara (2016) O. Palamara (ArgoNeuT),  _Proceedings, 10th International Workshop on Neutrino-Nucleus Interactions in the Few Region (NuInt15): Osaka, Japan, November 16-21, 2015_ , JPS Conf. Proc. 12, 010017 (2016).
  * 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] .
  * Rott  _et al._ (2017) C. Rott, S. In, J. Kumar, and D. Yaylali, JCAP 01, 016 (2017), arXiv:1609.04876 [hep-ph] .
  * Rott  _et al._ (2019) C. Rott, D. Jeong, J. Kumar, and D. Yaylali, JCAP 07, 006 (2019), arXiv:1903.04175 [astro-ph.HE] .
  * Serber (1947) R. Serber, Phys. Rev. 72, 1114 (1947).
  * Metropolis  _et al._ (1958) N. Metropolis, R. Bivins, M. Storm, J. 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. Bertsch, H. Kruse, and S. 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. Bertsch, W. Cassing, and U. Mosel, Phys. Rev. C 34, 2127 (1986).
  * Bertsch and Das Gupta (1988) G. Bertsch and S. Das Gupta, Phys. Rept. 160, 189 (1988).
  * Danielewicz and Bertsch (1991) P. Danielewicz and G. Bertsch, Nucl. Phys. A 533, 712 (1991).
  * Kadanoff and Baym (1962) L. Kadanoff and G. Baym,  _Quantum Statistical Mechanics_ (Benjamin, New York, 1962).
  * 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).
  * 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 .
  * 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] .
  * Mosel (2019) U. Mosel, J. Phys. G 46, 113001 (2019), arXiv:1904.11506 [hep-ex] .
  * Cugnon  _et al._ (1997) J. Cugnon, C. Volant, and S. Vuillier, Nucl. Phys. A 620, 475 (1997).
  * Duarte (2007) H. Duarte, Phys. Rev. C 75, 024611 (2007).
  * 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).
  * 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).
  * Glauber  _et al._ (1959) R. Glauber _et al._ , (1959).
  * Golubeva  _et al._ (1998) Y. Golubeva, L. Kondratyuk, A. Bianconi, S. Boffi, and M. Radici, Phys. Rev. C 57, 2618 (1998), arXiv:nucl-th/9712040 .
  * Uozumi  _et al._ (2012) Y. Uozumi, T. Yamada, S. Nogamine, and M. Nakano, Phys. Rev. C 86, 034610 (2012).
  * Boudard  _et al._ (2002) A. Boudard, J. Cugnon, S. Leray, and C. Volant, Phys. Rev. C66, 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 .
  * Battistoni  _et al._ (2013) G. Battistoni, F. Cerutti, A. Ferrari, J. Ranft, S. Roesler, and P. Sala, J. Phys. Conf. Ser. 408, 012051 (2013).
  * Battistoni  _et al._ (2015) G. Battistoni _et al._ , Annals Nucl. Energy 82, 10 (2015).
  * Golan  _et al._ (2012) T. Golan, C. Juszczak, and J. T. Sobczyk, Phys. Rev. C86, 015505 (2012), arXiv:1202.4197 [nucl-th] .
  * Andreopoulos  _et al._ (2010) C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Niewczas and Sobczyk (2019) K. Niewczas and J. T. Sobczyk, Phys. Rev. C 100, 015505 (2019), arXiv:1902.05618 [hep-ex] .
  * Carlson  _et al._ (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. Schmidt, and R. Wiringa, Rev. Mod. Phys. 87, 1067 (2015), arXiv:1412.3081 [nucl-th] .
  * Wiringa  _et al._ (1995) R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, Phys. Rev. C51, 38 (1995), arXiv:nucl-th/9408016 [nucl-th] .
  * Pieper (2008) S. C. Pieper, AIP Conference Proceedings 1011, 143 (2008), http://aip.scitation.org/doi/pdf/10.1063/1.2932280 .
  * Landau and Pomeranchuk (1953) L. D. Landau and I. Pomeranchuk, Dokl. Akad. Nauk Ser. Fiz. 92, 735 (1953).
  * Stodolsky (1975) L. Stodolsky, in  _Proceedings: International Colloquium on Multiparticle Reactions, 6th, Oxford, Eng., Jul 14-19, 1975_ (1975) p. 577.
  * Bogner  _et al._ (2010) S. Bogner, R. Furnstahl, and A. Schwenk, Prog. Part. Nucl. Phys. 65, 94 (2010), arXiv:0912.3688 [nucl-th] .
  * Barrett  _et al._ (2013) B. R. Barrett, P. Navratil, and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131 (2013).
  * Hagen  _et al._ (2014) G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. Dean, Rept. Prog. Phys. 77, 096302 (2014), arXiv:1312.7872 [nucl-th] .
  * Barbieri and Carbone (2017) C. Barbieri and A. Carbone, “Self-consistent Green’s function approaches,” (2017) pp. 571–644, arXiv:1611.03923 [nucl-th] .
  * Wiringa  _et al._ (2000) R. B. Wiringa, S. C. Pieper, J. Carlson, and V. Pandharipande, Phys. Rev. C 62, 014001 (2000), arXiv:nucl-th/0002022 .
  * Lovato  _et al._ (2013) A. Lovato, S. Gandolfi, R. Butler, J. Carlson, E. Lusk, S. C. Pieper, and R. Schiavilla, Phys. Rev. Lett. 111, 092501 (2013), arXiv:1305.6959 [nucl-th] .
  * Cruz-Torres  _et al._ (2019) R. Cruz-Torres, D. Lonardoni, R. Weiss, N. Barnea, D. Higinbotham, E. Piasetzky, A. Schmidt, L. Weinstein, R. Wiringa, and O. Hen, (2019), arXiv:1907.03658 [nucl-th] .
  * Sjostrand and van Zijl (1987) T. Sjostrand and M. van Zijl, Phys. Rev. D36, 2019 (1987).
  * (63) SAID database, http://gwdac.phys.gwu.edu/.
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Norbury (2008) J. Norbury, NASA Technical Paper 215116 (2008).
  * Dicello and Igo (1970) J. F. Dicello and G. Igo, Phys. Rev. C2, 488 (1970).
  * Bobchenko  _et al._ (1979) B. M. Bobchenko _et al._ , Sov. J. Nucl. Phys. 30, 805 (1979), [Yad. Fiz.30,1553(1979)].
  * Bauhoff (1986) W. Bauhoff, Atom. Data Nucl. Data Tabl. 35, 429 (1986).
  * Kox  _et al._ (1987) S. Kox _et al._ , Phys. Rev. C35, 1678 (1987).
  * Sihver  _et al._ (1993) L. Sihver, C. H. Tsao, R. Silberberg, T. Kanai, and A. F. Barghouty, Phys. Rev. C47, 1225 (1993).
  * Carlson (1996) R. F. Carlson, Atom. Data Nucl. Data Tabl. 63, 93 (1996).
  * Pandharipande and Pieper (1992) V. Pandharipande and S. C. Pieper, Phys. Rev. C 45, 791 (1992).
  * González-Jiménez  _et al._ (2019) R. González-Jiménez, A. Nikolakopoulos, N. Jachowicz, and J. Udías, Phys. Rev. C 100, 045501 (2019), arXiv:1904.10696 [nucl-th] .
  * O’Neill  _et al._ (1995) T. G. O’Neill _et al._ , Phys. Lett. B351, 87 (1995), arXiv:hep-ph/9408260 [hep-ph] .
  * Garrow  _et al._ (2002) K. Garrow _et al._ , Phys. Rev. C 66, 044613 (2002), arXiv:hep-ex/0109027 .
  * Garino  _et al._ (1992) G. Garino _et al._ , Phys. Rev. C 45, 780 (1992).
  * Radici  _et al._ (1994) M. Radici, S. Boffi, S. C. Pieper, and V. Pandharipande, Phys. Rev. C 50, 3010 (1994), arXiv:nucl-th/9408015 .
  * Benhar  _et al._ (1993) O. Benhar, V. Pandharipande, and S. C. Pieper, Rev. Mod. Phys. 65, 817 (1993).
  * Nikolaev  _et al._ (1993) N. N. Nikolaev, A. Szczurek, J. Speth, J. Wambach, B. Zakharov, and V. Zoller, Phys. Lett. B 317, 281 (1993), arXiv:nucl-th/9306007 .
  * Benhar (2013) O. Benhar, Phys. Rev. C 87, 024606 (2013), arXiv:1301.3357 [nucl-th] .
  * Benhar  _et al._ (1995) O. Benhar, A. Fabrocini, S. Fantoni, S. Pieper, V. Pandharipande, and I. Sick, Phys. Lett. B 359, 8 (1995).
  * Salcedo  _et al._ (1988) L. Salcedo, E. Oset, M. Vicente-Vacas, and C. Garcia-Recio, Nucl. Phys. A 484, 557 (1988).
  * Ashery  _et al._ (1981) D. Ashery, I. Navon, G. Azuelos, H. Walter, H. Pfeiffer, and F. Schleputz, Phys. Rev. C 23, 2173 (1981).
  * Ashery  _et al._ (1984) D. Ashery _et al._ , Phys. Rev. C 30, 946 (1984).
  * Jones  _et al._ (1993) M. Jones _et al._ , Phys. Rev. C 48, 2800 (1993).
  * Benhar  _et al._ (1994) O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A 579, 493 (1994).
  * Rocco  _et al._ (2016) N. Rocco, A. Lovato, and O. Benhar, Phys. Rev. Lett. 116, 192501 (2016), arXiv:1512.07426 [nucl-th] .
  * Benhar  _et al._ (2008) O. Benhar, D. day, and I. Sick, Rev. Mod. Phys. 80, 189 (2008), arXiv:nucl-ex/0603029 .
  * Cruz-Torres  _et al._ (2020) R. Cruz-Torres _et al._ (Jefferson Lab Hall A Tritium), Phys. Rev. Lett. 124, 212501 (2020), arXiv:2001.07230 [nucl-ex] .
  * Duer  _et al._ (2019) M. Duer _et al._ (CLAS), Phys. Rev. Lett. 122, 172502 (2019), arXiv:1810.05343 [nucl-ex] .

##  Appendix A Nucleon-nucleon interaction Models

The nucleon-nucleon cross sections used in our INC are obtained externally.
The elastic cross section is taken from the SAID database SAID database ,
obtained using GEANT4 Agostinelli  _et al._ (2003). The total cross section,
including inelastic contributions (e.g. pion production) is taken from Ref.
Norbury (2008). For reference, we reproduce those in Figs. 10. As it can be
seen, the cross sections from these databases are similar in the region where
the elastic contribution dominates, that is, below 300 MeV for $pp$
interactions and 500 MeV for $np$ interactions. The differences at high energy
can be attributed to inelastic contributions.

Figure 10: Comparison between the two interaction models used within the
intranuclear cascade model. The GEANT model is based on the SAID database SAID
database  and obtained using GEANT4 Agostinelli  _et al._ (2003) and the NASA
model is given from Norbury (2008).


# Towards Coalgebras in Stylometry

Joël A. Doat

###### Abstract

The syntactic behaviour of texts can highly vary depending on their contexts
(e.g. author, genre, etc.). From the standpoint of stylometry, it can be
helpful to objectively measure this behaviour. In this paper, we discuss how
coalgebras are used to formalise the notion of behaviour by embedding
syntactic features of a given text into probabilistic transition systems. By
introducing the behavioural distance, we are then able to quantitatively
measure differences between points in these systems and thus, comparing
features of different texts. Furthermore, the behavioural distance of points
can be approximated by a polynomial-time algorithm.

_K_ eywords Stylometry $\cdot$ Genre classification $\cdot$ Probabilistic
transition system $\cdot$ Behavioural distance $\cdot$ Coalgebra

##  1 Introduction

In stylometry, syntactic features of a text can be used as an indicator for
the author’s style. In this research, a new method to objectively measure this
behaviour is introduced. First, by embedding features into probabilistic
transition systems (PTSs) Larsen and Skou (1991), one can obtain a statistical
representation Sanderson and Günter (2006). These syntactic features include
for example co-occurrence of words or word classes, parts-of-speech
distribution, and grammatical structures Guillén-Nieto et al. (2008). Then,
the notion of behaviour in PTSs can be formalised by so-called probabilistic
bisimulations Larsen and Skou (1991). These are used to contrast the
probabilities of outgoing edges of two states in PTSs and return a value
describing the alikeness of these states, i.e. a lower value indicates greater
similarity.
When measuring the syntactical behaviour of a text with respect to a given
feature, the idea is to associate a state of the corresponding PTS as a
representant of the investigated text and compare it with the corresponding
state in the system of another text. Thus, for every feature, we obtain a
behavioural value for this text.
In order to calculate the value of probabilistic bisimulations, PTSs can be
translated into coalgebras Jacobs (2016) on metric spaces 111spaces with
distance measures Vink and Rutten (1999). Coalgebras can be seen as abstract
dynamic transition systems consisting of a state space (metric space) and a
transition function, describing one-step transitions at each state. The
advantage of this representation is that PTSs admit a so-called terminal
coalgebra Turi and Rutten (1998) which intuitively describes a transition
system realising all possible state behaviours while being a metric space
itself again. In other words, when investigating texts, the terminal coalgebra
can be seen as an abstract ontological space containing information (about the
behaviour with respect to a given feature) of every possible text together
with a distance measure. Thus, we can measure the so-called behavioural
distances Breugel and Worrell (2006) between different behavioural
realisations on the terminal coalgebra. Furthermore, by definition, for every
state, there is then a unique map that associates it with its behavioural
realisation. The map can be approximated by a polynomial-time algorithm and
thus, the behavioural value of the investigated state. More specifically, this
algorithm solves a transportation problem Chvátal (1983) that can be reduced
to a minimum cost network flow Orlin (1988) and minimum cost circulation
Tardos (1985) problem.
At last, this algorithm will be tested to obtain the first results. For our
purpose, chosen texts will be analysed in contrast to the genre categories of
the Brown corpus.

##  2 Method

###  2.1 Probabilistic Transition System

Probabilistic transition systems can be used as a method to represent
statistical data with respect to the influence of its environment. The general
idea is that two information states (e.g. lexical items in a text) can be
related via a transitioning probability and a provoking action. Recall the
formal definition.

######  Definition 1.

Larsen and Skou (1991) A probabilistic transition system (PTS) is a triple
$(S,A,\pi)$ with a non-empty set of states $S$, a set of actions $A$ and a
transition function $\pi:S\times A\times S\to[0,1]$ such that for all $s\in S$
and $a\in A$, $\sum_{s^{\prime}\in S}\pi(s,a,s^{\prime})\in[0,1]$.

Note that PTSs without additional labels on the edges (unlabelled case) can be
defined by a set of actions only containing an arbitrary element (i.e.
$A:=\\{\ast\\}$). This allows to make no distinction between edges with
respect to actions. For notational simplicity, we only consider the unlabelled
case and denote PTSs as a tuple $(S,\pi)$.

Now, consider the following texts.

By labelling the edges of the corresponding word adjacency networks (WAN) with
the transition probabilities, we obtain the following diagram. Note that the
blue and and red arrows belong to the left and right text, respectively.

Figure 1: Probabilistic Transition System for Text A and B

As one can see in Fig.1, using probabilistic transition systems (PTSs), the
lexical unit ’can’ emits a different behaviour depending on the context which
can be chosen in terms of linguistic features that can be embedded in PTSs.
The resulting system represents a ’probabilistic reading flow’ with respect to
the given features. The question might arise by ’how much’ the behaviour
differs for a specific word.

For the unlabelled case (no additional actions), such models are already in
use for stylometric analysis, e.g. sequence of characters Khmelev and Tweedie
(2008), sequence of words Sanderson and Günter (2006) and sequence of function
words Eisen et al. (2018). Usually, this is done by averaging values of a
transition matrix given by a text and comparing different texts with these
results. Since these systems are used to represent the syntactic style of a
text, the corresponding behaviour can be seen as a core property of this style
representation.

###  2.2 Probabilistic Bisimulation

In this context, the qualitative equivalence notion to check if emitted
behaviour between states of a PTS coincide is called probabilistic
bismulation.

######  Definition 2.

Larsen and Skou (1991) For a set of states $S$ and transition function $\pi$
define a probabilistic transtion system $(S,\pi)$. An equivalence relation
$\mathcal{R}$ on the set of states $S$ is a probabilistic bisimulation if
$s\mathcal{R}s^{\prime}$ implies $\sum_{\bar{s}\in
E}\pi(s,\bar{s})=\sum_{\bar{s}\in E}\pi(s^{\prime},\bar{s})$ for all
$\mathcal{R}$-equivalence classes $E$. States $s$ and $s^{\prime}$ are
probabilistic bisimilar if $s\mathcal{R}s^{\prime}$ for some probabilistic
bisimulation $\mathcal{R}$.

Such a bisimulation is intended to be a witness for equivalent behaviour
between states, i.e. two states are equivalent if they have the same
transitioning probability to every equivalence class of the bisimulation.

Consider the follwoing PTS.

Figure 2: Probabilistic Transition System t compare behaviour of $w$ and $c$

In Fig.2, the states $w$ and $c$ emit the same iff $\varepsilon=1/6$. Thus, in
this simple case, we would like to define the difference between $w$ and $c$
as $1/6$ if $\varepsilon=0$. Motivated by this definition, Breugel and Worrell
(2006) use a pseudometric222A pseudometric is like a usual metric with the
difference that two distinct elements can have distance 0. on states of PTSs
as a method to witness a quantitative version of this behavioural equivalence,
i.e. a smaller distance signifies greater similarity. The idea is to represent
these PTSs as coalgebras Jacobs (2016) on metric spaces333More precisely, PTSs
can be seen as coalgebras on the category of $1$-bounded complete ultrametric
spaces and non-expansive functions Vink and Rutten (1999).. The advantage of
this representation is that PTSs admit a so-called terminal coalgebra444See
Rutten and Turi’s ultrametric terminal coalgebra theorem Turi and Rutten
(1998) which intuitively describes a transition system realising all possible
state behaviours, while being a metric space itself again. By definition of
terminal coalgebras, for every state there is then a unique map that
associates it with its behavioural realisation. The distance between such
realisations is called behavioural distance.

###  2.3 Coalgebra

In this subsection, I introduce coalgebras and the behavioural distances in a
rather intuitive way to focus on how this idea interacts with stylometry
instead of introducing all necessary terms from category theory.

######  Definition 3.

Jacobs (2016) For an (set) endofunctor $T$, object (set) $C$, and morphism
$\gamma:C\to TC$ the tuple $(C,\gamma)$ is called a $T$-coalgebra.

Coalgebras can be seen as abstract state-based transition systems. It is a
more general notion than PTSs because it does not restrict to a set of states
with transitioning probabilities between them. Such a $T$-coalgebra
$(C,\gamma)$ consists of three parameters that describe it.

  1. •

The object (set) $C$ is called the carrier of a coalgebra and denotes the
states of the transition system. In the case of PTSs, $C$ is a set of states.

  2. •

The type $T$ describes all observable behaviours after one step. Applying this
endofunctor transforms the object of a type (e.g. set) into an object of the
same type. For PTSs, $T$ is the so-called distribution functor which maps a
set of states $C$ to a set $TC$ of all possible probabilistic distributions on
the same states. $TC$ is also referred to as structured successor.

  3. •

$\gamma$ denotes the actual transition function of the system, i.e. mapping
one element of the carrier to the behaviour it is emitting. In PTSs, this is a
function mapping a state to its corresponding probabilistic distribution.

Since the set $TC$ contains all possible observable behaviours of states in
$C$ after one step, we can apply $T$ again and obtain with $T^{2}C$ the
observable behaviours after two steps. This process can be continued
iteratively to describe the behaviour after arbitrarily many steps.
As mentioned in the last subsection, the existence of a terminal coalgebra
guarantees that the information about the behaviour of $(C,\gamma)$ collapses
through the iteration process into a single point inside the terminal
coalgebra itself. In other words, the terminal coalgebra realises every
possible $n$-step behaviour with respect to $T$ and thus, there is one state
that emits the same behaviour as a state in the coalgebra we started with.
Therefore, states in the terminal coalgebra can be seen as behavioural
realisations of transitions systems.

Figure 3: Observable behaviour collapses to one point in the terminal
coalgebra

As PTSs are special instances of such abstract transition systems, we can
embed the statistical representation of texts (see Fig.1) into a coalgebra. In
Fig.3, ($C$,$\gamma$) and $(C^{\prime},\gamma^{\prime})$ constitute now two
texts as coalgebras. For a chosen starting point in the PTSs the observable
behaviours of these texts with respect to the representation collapse to point
on the terminal coalgebra. To measure the behavioural distance which
correspond to a stylistic distance in our context, we need a metric on the
terminal coalgebra. Breugel and Worrell (2006) showed that PTSs can be
equivalently represented as transition systems on pseudometric spaces such
that also the distance of probabilistic distributions in the structured
successors can be measured. The advantage of this is that terminal coalgebra
also becomes a pseudeometric space and thus, behavioural/stylistic distances
between texts exist. In other words, when investigating texts, this terminal
coalgebra is a naive approach to build an ontology containing
behavioural/stylistic information (with respect to a given syntactic features)
of every possible text together with a distance measure.

###  2.4 Algorithm

On the basis of the idea discussed in Section 2.3, Breugel and Worrell (2006)
extracted a polynomial algorithm that approximates the value of behavioural
distances. The algorithm solves a linear programming problem.

Figure 4: Network flow model of the transportation problem

In a PTS with $N$-many states, the outgoing probabilities of two states
$p_{1},\dots,p_{N}$ and $q_{1},\dots,q_{N}$, respectively, are compared by
solving a transportation problem. In Fig.4, one can see the corresponding
network flow model. The idea is to find the minimal transportation costs
$c_{i,j},i,j=1,\dots,N$ from the source $q_{i}$ to the destination $p_{j}$
where the outgoing probabilities label the states with their supply (source)
and demand (destination). To encapsulate also more than one transition step
and thus, increase the accuracy of the result, this process can be iterated
arbitrarily many times555One can imagine the iteration process in each step as
an additional comparison along the observable $n$-step behaviour as shown in
Fig.3.. Furthermore, thanks to Proposition 12 in Breugel and Worrell (2006),
the accuracy $\alpha$ can be chosen freely since for some constant
$c\in(0,1)$666The value of this constant is prescribed by the proposition. the
difference between the value of the behavioural distance at step $\lceil
log_{c}(\alpha/2)\rceil$ and the actual distance $d$ is smaller or equal to
$\alpha$. Note that the calculated result is then $c\cdot d$.

// Step 1: Initialisation

foreach  _$1\leq k,l\leq N$_ do

$d_{kl}\leftarrow 0$

end foreach

$c_{00}\leftarrow 0$

foreach  _$1\leq i,j\leq N$_ do

$c_{i0}\leftarrow 1$

$c_{0j}\leftarrow 1$

end foreach

// Step 2: Main loop

repeat  _$\lceil log_{c}(\alpha/2)\rceil$_ times

foreach  _$1\leq i,j\leq N$_ do

$c_{ij}\leftarrow c\cdot d_{ij}$

end foreach

foreach  _$1\leq k,l\leq N$_ do

/* min_value defined below */

$d_{kl}\leftarrow\text{min\\_value}(k,l)$

end foreach

end

Figure 5: Pseudocode of the algorithm

Fig.5 describes the algorithm in pseudocode. For a PTS $(S,\pi)$ with states
$S=\\{s_{1},\dots,s_{N}\\}$, the operation min_value$(k,l)$ is defined in each
iteration step by the minimal value of

$\displaystyle\sum_{i,j=0}^{N}c_{ij}\lambda_{ij}$

with constraints

$\displaystyle\sum_{i,j=0}^{N}\lambda_{ij}=\pi(s_{k},s_{j}),$ $\displaystyle
0\leq j\leq N,$ $\displaystyle\sum_{i,j=0}^{N}\lambda_{ij}=\pi(s_{l},s_{i}),$
$\displaystyle 0\leq i\leq N,$ $\displaystyle 0\leq\lambda_{ij},\,0\leq
i,j\leq N,$

where

$\displaystyle\pi(s_{i},s_{0}):=1-\sum_{j=1}^{N}\pi(s_{i},s_{j}).$

The resulting matrix $d$ is the distance matrix, i.e. the distance between
states $s_{k}$ and $s_{l}$ is the value $d_{kl}$.

##  3 Application

###  3.1 Idea

To be able to use the algorithm on two different texts, we need to combine the
corresponding PTSs into the same system. As we saw in Fig.2, the idea of
behavioural distance is to compare transition probabilities to the same
equivalence classes of states. This means that the easiest way to compare
texts is to connect them at some chosen end state. In other words, all
transition walks of both PTSs can be extended to a walk whose final node
corresponds to the designated end state and the comparison of states
materialise in terms of their relation to it.
Secondly, we need states that represent the whole text in each PTS such that
they can be juxtaposed to measure a stylistic distance.

Figure 6: How to compare two texts

Fig.6 demonstrates an example on how to apply this method for the purpose of
stylometry. Given texts $T_{1}$ and $T_{2}$, one could measure the difference
between the distribution of Parts-of-Speech (POS), i.e. a PTS with POS as
states and the transitioning probabilities given by their co-occurrences in
the text. In order to combine both PTSs into one, the punctuation marks at the
end of the sentences are represented by a single state which coincides for the
two texts. Each state representing a POS has then a transition to this end
state labelled by the probability that a sentence ends after given POS. This
allows a clear comparison by using the end state as the only state in the
equivalence class accessible by both systems. Next, we add a starting state in
each PTS that has transitions to each state representing a POS that occur at
the beginning of a sentence and label these transitions with the probability
that a sentence starts with the given POS. The two starting states can then be
compared in terms of the behavioural distance. The resulting value represents
the syntactical difference for the POS distributions.
Further linguistic features/co-occurrence phenomenons that can be represented
as such a transition system include: lexical items, letters, collocations of
specific POS (e.g. nouns after modal verbs), function words, and grammar trees
of sentences. In the next section, I used multiple features to obtain a more
precise result. More specifically, after the same principle as in Fig.6, I
also used the letter distribution (states represent letters in the observation
instead of POS) and the grammar trees (states represent grammatical
constituents in the hierarchy). The corresponding starting states connect to
the first letters in sentences and the first separation into phrasal
constituents in the grammatical hierarchy of sentences, respectively. The end
state is given by punctuation marks and the lowest level in the grammar
trees777For both PTSs that represent the statistical relation of grammar trees
in the texts, the last level in the hierarchy corresponds to the actual
lexical items that are irrelevant for the grammatical observation. Thus, all
levels consisting of lexical items collapse to the end state of the combined
system., respectively.

###  3.2 Results

To demonstrate the idea, I applied the algorithm in the context of genre
classification. For this purpose, we take 4 texts and compare each of them
with texts that represent genres as discussed in the last subsection (see
Fig.6). The genres are news, religion, government, and belles-lettres from the
Brown corpus. The texts are the news about coconuts (Reuters corpus), the
bible (Gutenberg corpus), Donald Trump’s inauguration speech (Inaugural
Address corpus), and Lewis Carroll’s Alice’s Adventures in Wonderland
(Gutenberg Corpus). For all comparisons we investigate the syntactic features
Parts-of-Speech, grammar trees, and letters. The three resulting values can be
seen as a vector whose distance to $0$ can be calculated with the euclidean
distance in order to compare it with other categories. The results can be
found in Fig.7-10.

|  news |  religion |  gov. |  belles
---|---|---|---|---
POS |  0.2778 |  0.3899 |  0.3689 |  0.41
Grammar |  0.322 |  0.4111 |  0.3556 |  0.4444
Letters |  0.083 |  0.0967 |  0.0877 |  0.1222
Euclid |  0.4333 |  0.5748 |  0.5199 |  0.6166
Figure 7: Results for coconut news |  news |  religion |  gov. |  belles
---|---|---|---|---
POS |  0.3661 |  0.2795 |  0.2633 |  0.314
Grammar |  0.4333 |  0.2919 |  0.3201 |  0.4005
Letters |  0.1007 |  0.0501 |  0.0199 |  0.0885
Euclid |  0.5761 |  0.4072 |  0.415 |  0.5166
Figure 8: Results for the bible |  news |  religion |  gov. |  belles
---|---|---|---|---
POS |  0.2519 |  0.2238 |  0.2115 |  0.2903
Grammar |  0.2581 |  0.273 |  0.2667 |  0.3856
Letters |  0.0447 |  0.0329 |  0.0454 |  0.0548
Euclid |  0.359 |  0.3545 |  0.3434 |  0.4858
Figure 9: Results for Donald Trump’s inauguration speech |  news |  religion |  gov. |  belles
---|---|---|---|---
POS |  0.2868 |  0.2175 |  0.2284 |  0.2273
Grammar |  0.3781 |  0.3754 |  0.3502 |  0.3389
Letters |  0.0732 |  0.0459 |  0.0177 |  0.0638
Euclid |  0.4802 |  0.4363 |  0.4185 |  0.413
Figure 10: Results for Alice’s Adventures in Wonderland

As one can see in these results, the distance between a text and its correct
category is the lowest. Thus, for this small experiment the classification
task worked. In the context of this demonstration, it is clear that the method
of behavioural distance gives an interesting research direction for
stylometry. Nevertheless, further experiments are still required to check if
this method is reliable since the difference between some distances is rather
small (see euclidean distances to government and belles-lettres in Fig.10).

##  4 Conclusion

In this paper, I introduced a new tool for objective measurement of text
behaviour. By embedding syntactic features of texts into PTSs and using the
method of Breugel and Worrell (2006), we have shown how to calculate the
distance of behaviour realisations. Furthermore, a polynomial-time algorithm
exists to approximate this distance and according to the chosen syntactic
feature, it delivers a vector-like representation of the investigated text. In
a small experiment we have seen a small task of genre classification which has
been solved correctly with this method.
There are several aspects that can be improved and are subject to future work.
So far, I only used the features regarding Parts-of-Speech, grammar trees, and
letters. As already mentioned, further features that can be represented as
PTSs include lexical items, collocations of specific POS (e.g. nouns after
modal verbs), and function words.
In Fig.6, we saw one idea of how two PTSs can be combined. Since this
combination is variable, it remains to be investigated if there are better
ways to achieve this. Another direction one could possibly go is to embed a
PTS into a canonical construction instead, i.e. extending every PTS uniformly
into a bigger system such that every text obtain a behavioural distance value
with respect to the same construction without a second text.
Two properties of PTSs we did not make use of in this paper are actions and
subprobabilities.

Figure 11: Labels and subprobabilities

Fig.11 demonstrates a simple example of such a case. The colors green and
orange denote two different actions on $s_{0}$ that have distinct probability
distributions. According to Breugel and Worrell (2006) the algorithm can also
be extended for actions. Furthermore, the sum of all probabilities along the
orange action is smaller than $1$. This can be interpreted as an additional
probability that the system refuses to interact with the environment at
$s_{0}$. By adding states representing punctuation marks like comma, colon,
and semicolon, one can include subprobability distributions at these states to
denote an interruption in the reading flow.
An open question not covered in this paper is the possible use of this method
in other contexts like for example pragmatics or register analysis. A state
representing a concrete lexical item can be introduced and the corresponding
PTS describes a co-occurrence phenomenon (e.g. parts-of-speech occurring after
the modal verb "can"). This can be useful measure the behavioural distance of
this item in different contexts.

## References

  * Breugel and Worrell (2006) Breugel and Worrell.  Approximating and computing behavioural distances in probabilistic transition systems.  _Theoretical Computer Science_ , 360(1-3):373–385, 2006.
  * Chvátal (1983) Chvátal.  _Linear Programming_.  W.H. Freeman and Company, New York, San Francisco, 1983.
  * Eisen et al. (2018) Eisen, Ribeiro, Segarra, and Egan.  Stylometric analysis of early modern period english plays.  _Digital Scholarship in the Humanities_ , 33(3):500–528, 2018.
  * Guillén-Nieto et al. (2008) Guillén-Nieto, Vargas-Sierra, Juan, Martinez-Barco, and Suárez Cueto.  Exploring state-of-the-art software for forensic authorship identification.  _Journal of English Studies_ , 8(1), 2008.
  * Jacobs (2016) Jacobs.  _Introduction to Coalgebra: Towards Mathematics of States and Observation_.  Cambridge University Press, Cambridge, UK, 2016.
  * Khmelev and Tweedie (2008) Khmelev and Tweedie.  Using markov chains for identification of writers.  _Literary and Linguistic Computing_ , 16(4):299–307, 2008.
  * Larsen and Skou (1991) Larsen and Skou.  Bisimulation through probabilistic testing.  _Information and Computation_ , 94(1):1–28, 1991\.
  * Orlin (1988) Orlin.  A faster strongly polynomial minimum cost flow algorithm.  In _Proceedings of the 20th Annual ACM Symposium on Theory of Computing_ , pages 377––387, Chicago, USA, 1988.
  * Sanderson and Günter (2006) Sanderson and Günter.  Short text authorship attribution via sequence kernels, markov chains and author unmasking: An investigation.  In _Proceedings of the 2006 Conference on Empirical Methods in Natural Language Processing_ , pages 482–491, Sydney, Australia, 2006. Association for Computational Linguistics.
  * Tardos (1985) Tardos.  A strongly polynomial minimum cost circulation algorithm.  _Combinatorica_ , 5(3):247––255, 1985.
  * Turi and Rutten (1998) Turi and Rutten.  On the foundations of final semantics: non-standard sets, metric spaces, partial orders.  _Mathematical Structures in Computer Science_ , 8(5):481––540, 1998.
  * Vink and Rutten (1999) De Vink and Rutten.  Bisimulation for probabilistic transition systems: a coalgebraic approach.  _Theoretical Computer Science_ , 221(1-2):271–293, 1999.