Slow neutron production as a probe of hadron formation in high-energy $\gamma^*A$ reactions

TL;DR

This paper investigates low-energy neutron production in high-energy photon and muon-induced nuclear reactions, using hybrid models to understand hadron formation and final state interactions, with implications for future collider experiments.

Contribution

It provides new hybrid dynamical+statistical calculations of neutron production in $ u$- and photon-nucleus collisions, linking experimental data to hadron formation models and small-$x$ dynamics.

Findings

E665 data shows strong suppression of final state interactions for hadrons above 1 GeV/c.

Neutron multiplicity and $p_t$ spectra predictions for LHC and RHIC energies.

Ultraperipheral heavy-ion collisions can probe small-$x$ dynamics and photon wave function components.

Abstract

Deep Inelastic Scattering (DIS) experiments at the planned Electron-Ion Collider will be affected by details of the hadron formation inside the nuclear volume. Besides semi-inclusive particle production experiments decays of the target nucleus via emission of neutrons provide an additional opportunity to probe this domain. This paper reports on the hybrid dynamical+statistical calculations of low-energy neutron production in muon- and virtual photon-induced collisions with nuclei. We confirm the conclusion that the E665 data on neutron production in $\mu^-$ + Pb DIS at 470 GeV indicate a strong suppression of the final state interaction for hadrons with momenta above $\sim 1$ GeV/c. Ultraperipheral heavy-ion collisions at the Large Hadron Collider (LHC) and the Relativistic Heavy Ion Collider (RHIC) can be used to test this suppression. The calculations of the neutron multiplicity distributions and $p_t$-spectra in photon - nucleus collisions at the energies accessible at the LHC and RHIC are presented for several models of hadron formation. We argue that studies of neutron production in ultraperipheral heavy ion collisions open a new window on the small-$x$ dynamics and hadron component of the photon wave function.

Tables

Table 1. Table 1: Parameters of UPC Au+Au at RHIC and Pb+Pb at LHC.

	$\sqrt{s_{NN}}$ (TeV)	${\gamma }_L$	$k^{\max }$ (TeV/c)	$W$ (GeV)
RHIC	0.2	106	0.642	34.7
LHC	5.5	2931	477	946

Equations

$$d{P}_i=\frac{{\rho }_i(\mathbf{r})d^{3}r}{A},$$

$$\frac{d\sigma }{d\Omega d{E}^{\prime }}=\Gamma [{\sigma }_T(W^2,Q^2)+ϵ{\sigma }_L(W^2,Q^2)],$$

$$\Gamma =\frac{\alpha E^{\prime }(W^2-m_N^2)}{(2\pi {)}^2{Q}^2{m}_{N}E(1-ϵ)}$$

$$ϵ$$

$$\frac{d\sigma }{dyd{Q}^2}=\frac{\pi }{E^{\prime }}\frac{d\sigma }{d\Omega d{E}^{\prime }},$$

$$l_{\mathrm{form}}=a{p}_h,$$

$${\sigma }_{\mathrm{eff}}(t)/{\sigma }_0=X_0+(1-X_0)\frac{t-t_{\mathrm{prod}}}{t_{\mathrm{form}}-t_{\mathrm{prod}}},$$

$${\sigma }_{\mathrm{eff}}/{\sigma }_0=\Theta (p_{\mathrm{cut}}-p_h),$$

$$A_{\mathrm{res}}$$

$$Z_{\mathrm{res}}$$

$$E_{\mathrm{res}}^{\ast }$$

$${\mathbf{p}}_{\mathrm{res}}$$

$$n_{\mathrm{wound}}=1/R_A\sim 1.25,$$

$$x_F=\frac{E-p^z}{(E_A-p_A^z)/A},$$

$$x_F=\frac{E-p^z}{m_N}.$$

Full text

1 Introduction

The planned Electron-Ion Collider (EIC) will give us access to the quark- and gluon-substructure of bound nucleons in energy regimes that extend the earlier experiments at Hadron-Electron Ring Accelerator (HERA) and Thomas Jefferson National Accelerator Facility (JLab) into a new domain Accardi et al. (2016). Experiments at the EIC will use the nucleus both as a target and a laboratory for quantum chromodynamics (QCD) where, e.g., the quark hadronization can be studied. In all these experiments one tacitly assumes an approximate factorization of the initial DIS process and the final state interactions (FSI) of the particles produced. All experimental observables will necessarily be affected by both. To gain any knowledge of the microscopic quark-gluon structure of bound nucleons it is, therefore, mandatory to have the FSI under quantitative control.

In the DIS processes, the production of a leading hadron involves a multistep fragmentation of a leading quark. The attenuation of leading hadrons has been measured by the European Muon Collaboration (EMC) Ashman et al. (1991) and HERMES Collaboration Airapetian et al. (2001). In the models of refs. Accardi et al. (2003); Akopov et al. (2002); Gallmeister and Falter (2005); Gallmeister and Mosel (2008), which are based on the string fragmentation, the attenuation is explained by the pre-hadronic absorption on surrounding nucleons. In the model of ref. Kopeliovich et al. (2004), both the gluon bremsstrahlung and the pre-hadronic absorption are included to describe the experimental data. In the models of refs. Wang and Wang (2002); Arleo (2003), the same data are described by modifications of the quark fragmentation function due to the medium-induced gluon radiation.

So far, the only theoretical analysis of both EMC Ashman et al. (1991) and HERMES Airapetian et al. (2001) experiments is that of Gallmeister et al. in ref. Gallmeister and Mosel (2008). These authors used a transport model to follow the initial fragmentation all the way to the asymptotic, final state of all ejected particles. This is necessary since experimental acceptances can overshadow the underlying microscopic mechanism as was illustrated by Falter et al. in ref. Falter et al. (2004).

At high energies and hard scattering processes color transparency, which reduces the FSI of initially produced fast particles, is expected to be present and can be studied in semi-inclusive lepton-induced particle production spectra. An alternative is the study of hadron production in the nuclear fragmentation region. If the produced fast particles interact with the target nucleus only weakly the excitation of the target remnant should be small. It is, therefore, promising to study the production of low-energy hadrons emitted from the target remnant. Such studies were performed so far only by the E665 Collaboration at Fermilab. In ref. Adams et al. (1995a), the production of low energy ($E<10$ MeV) neutrons in $\mu^{-}$+Pb DIS at 470 GeV have been studied. The first theoretical analysis of these data has been performed in ref. Strikman et al. (1999). These authors have demonstrated that the energy dependence of the neutron spectrum can be reproduced for a range of assumptions about the formation time. However, in order to reproduce the total neutron multiplicity of about five neutrons per event, it was necessary to assume that only nucleons produced in DIS with momenta $\leq$ 1 GeV/c could interact. This surprisingly low value is well below the momentum range where the formation time and color transparency effects become important. This hints at a new dynamics for the hadron formation in the nucleus fragmentation region.

In another E665 experiment Adams et al. (1995b) on $\mu^{-}$+Xe DIS at 490 GeV, the multiplicities of so-called gray tracks, ascribed to protons in the momentum range 200-600 MeV/c, have been measured. The theoretical analysis of this experiment in ref. Ciofi degli Atti and Kopeliovich (2005) has shown that the data on the $Q^{2}$- and $x$-dependence of the average number of gray tracks can be described when one assumes the color transparency only for the leading hadrons, in line with the assumptions of ref. Strikman et al. (1999). In this model only slow recoiling nucleons from elementary lepton-nucleon DIS re-interact with the nucleus.

In this paper, we suggest a way to study the formation of hadrons in the nucleus fragmentation region at the collider energies before the advent of the EIC. The proposed strategy is based on the experience gained from the studies of the ultraperipheral collisions (UPCs) of heavy ions at the LHC and direct photon interactions at HERA. In analogy with DIS, one can select a hard interaction of a quasireal photon coherently emitted by one of the nuclei with a parton of the other nucleus by triggering on the production of two high $p_{T}$ jets which together carry practically all energy of the photon. They have large invariant mass, are mostly produced at central rapidities and can be easily detected by the CMS and Atlas detectors Angerami (2018). The hard photon interaction thus initiates a cascade of particles. If the nucleus is more transparent for the energetic products, counting slow neutrons with energies up to about 10 MeV then uses the nucleus as a *”micro calorimeter” *

In the present study we simulate the initial $\gamma^{*}N$ event on the bound nucleon via the PYTHIA model 333We will mostly discuss the PYTHIA results in terms of photon-gluon interactions which are predominant in the studied region of Bjorken $x\sim 10^{-2}-10^{-1}$, although also a moderate admixture of photon-quark interactions can be present in PYTHIA events as well.. The slow neutrons are mostly produced in the interaction of the struck nucleon debris with surrounding nucleons and we keep track of these interactions by means of quantum-kinetic transport theory as encoded in the Giessen Boltzmann-Uehling-Uhlenbeck (GiBUU) model Buss et al. (2012). For the low-energy particles we supplement this description with the statistical multifragmentation model (SMM) Botvina et al. (1987); Bondorf et al. (1995) that describes the evaporation of low-energy neutrons by an excited nuclear residue.

The paper is organized as follows. Section 2 gives an overview of the theoretical framework containing the $\gamma^{*}N$ event simulation by the PYTHIA model included in the GiBUU package, the propagation of produced particles through the nucleus which leads to the energy deposition in the form of hole excitations in the nuclear residue, and the statistical deexcitation of the nuclear residue. Numerical results are presented in sec. 3. First, we benchmark the model against experimental data on slow neutron production in $\mu^{-}+A$ collisions at 470 GeV. Next, we perform calculations for the LHC and RHIC kinematics of ultraperipheral heavy-ion collisions in which one of the nuclei serves as a source of hard quasireal photons interacting with another nucleus. We show that the multiplicity distribution and $p_{t}$-spectrum of soft neutrons are directly related to pre-hadronic interactions in the nucleus. In sec. 4, we discuss directions for further studies, in particular, the expectations for enhanced neutron production in the resolved photon kinematics. Finally, in sec. 5 we summarize our results and draw some conclusions.

Appendix A contains the numerical results of simulations of proton-nucleus reactions at $p_{\rm lab}=1-2$ GeV/c as a test of the interface between the dynamical and statistical parts of our hybrid theoretical approach. We determine the time window when the nuclear residue reaches thermal equilibrium and thus the SMM becomes applicable. The energy spectra of neutrons from $pA$ collisions are well described by our hybrid calculations in a quite broad range of neutron energy, up to 160 MeV. Remarkably, the neutron yields below 10 MeV are governed by statistical evaporation and agree very well with experiment.

2 Theoretical framework

We are addressing slow particle production in the interaction of a highly-energetic (virtual) photon with a nuclear target. The physical picture of this process is quite complicated and, in-fact, is not fully covered by any existing theoretical approach. Hence, largely guided by experience gained from previous studies of the particle- and nucleus-nucleus interactions Barashenkov et al. (1974); Bondorf et al. (1995); Botvina et al. (1995); Strikman et al. (1999); Buss et al. (2012) we rely on the hybrid description combining three steps: (1) hard interaction of the incoming $\gamma^{*}$ with a bound nucleon, (2) propagation of produced particles in the target nucleus, (3) statistical evaporation of slow neutrons from a nuclear residue. For each step we use a dedicated and well tested in its energy range theoretical model. Below in this section, we describe all three model steps and – in some more detail – the interfaces between models, i.e. the hadron formation process and the determination of the parameters of the nuclear residue.

2.1 Hard interaction

In the case of high energy lepton-nucleus interaction the struck nucleon is chosen randomly neglecting nuclear shadowing effects, see discussions below in secs. 3.1 and 3.2. This corresponds to the probability distribution of the interaction with protons and neutrons

[TABLE]

where $\rho_{i}(\boldsymbol{r})~{}(i=p,n)$ are the densities of protons and neutrons. The collision of virtual photon with the struck nucleon is simulated via PYTHIA version 6.4 Sjostrand et al. (2006).

The kinematics of the lepton scattering depends on experimental conditions. In the case of DIS the elementary lepton-nucleon cross section can be written as follows (cf. Christy and Bosted (2010)):

[TABLE]

where $\Omega$ and $E^{\prime}$ are, respectively, the solid angle and the energy of the scattered lepton in the laboratory frame,

[TABLE]

is the flux of virtual photons, $\sigma_{T}(\sigma_{L})$ is the photoabsorption cross section for purely transverse (longitudinal) photons, and

[TABLE]

is the relative flux of longitudinal virtual photons, $y=\nu/E$ is the fraction of the beam energy $E$ carried off by a virtual photon. The second equality of Eq.(4), which is actually implemented in GiBUU, is satisfied in the limit of the zero lepton mass exactly.

In actual simulations, Eq.(2) has been applied in the rest frame of the struck nucleon (which has a finite momentum due to the Fermi motion) by using the relation

[TABLE]

where, in the Lorentz-invariant form, $y=pq/pk$ with $p,k$, and $q$ being the four-momenta of the struck nucleon, initial lepton, and the virtual photon, respectively. Note that in the present GiBUU calculations the reduced cross section $\sigma^{*}=\sigma_{T}+\epsilon\sigma_{L}$ is extracted from PYTHIA with default settings (MSTP(14)=30), i.e. it includes the mixture of all possible processes (point-like photon, VMD etc.).

2.2 Hadron formation

Particles produced in a hard interaction are not immediately formed and thus interact with the target nucleons not like ordinary hadrons with a fully developed internal structure. It is well accepted that the interaction of such pre-hadrons with target nucleons is significantly reduced. A pre-hadron converts to the ordinary hadron within the distance called formation length or, equivalently, within the time interval called formation time.

In ref. Gallmeister and Mosel (2008), the consistent description of the HERMES and EMC data has been reached within the GiBUU transport model with the pre-hadron production and formation times of the yo-yo picture. Following this picture, the difference between the formation and production proper times for a hadron with mass $m_{h}$ should be of the order of $m_{h}/2\kappa$ (1/4 of the $q\bar{q}$ string oscillation period Andersson et al. (1983)) where $\kappa=1$ GeV/fm is the string tension.

This estimate can be compared with that of the color transparency model for the hadron attenuation in exclusive processes like $eA\to e+\mbox{meson}+A^{*}$ at large enough $Q^{2}$ and momenta of the produced meson, see Dutta et al. (2013) for review and references. In this case, the color-neutral $q\bar{q}$ dipole of the transverse size $r_{\perp}\sim 1/Q$ is formed and interacts with the nucleon with the cross section $\propto r_{\perp}^{2}$ as follows from perturbative QCD. Such a “point-like” configuration is not static. However, it can be decomposed in the hadronic basis of the states with the same momentum but different energies. Therefore, the decomposition looses its coherence on the time scale $\sim 2p_{h}/\Delta M^{2}$ where $\Delta M^{2}$ is of the order of the difference between the squared masses of the particle to be formed and its nearest radially excited state. This time scale can be regarded as a formation time of the hadron, or the expansion time of the point-like configuration to the normal hadronic size. It is expected from the analyses of pion electroproduction at JLab Larson et al. (2006); Larionov et al. (2016), that $\Delta M^{2}\simeq 0.7$ GeV2. This corresponds to almost the same value of the hadron formation proper time as predicted by the yo-yo picture.

Thus, the hadron formation length is similar in DIS and exclusive processes and can be written as

[TABLE]

where $p_{h}$ is a hadron momentum, and $a=(0.4-0.6)$ fm c/GeV is a constant factor. Eq.(6) provides a guideline for our present-day understanding of the hadron formation length. This equation suggests that the hadron will be formed inside a heavy nucleus, i.e. the formation length is less than a nuclear radius, if p_{h}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}10 GeV/c, while at higher momenta the pronounced transparency should be observed.

To describe the hadron formation effects we applied the three different prescriptions:

(i) Default formation procedure of the GiBUU model Gallmeister and Falter (2005); Gallmeister and Mosel (2008); Buss et al. (2012) based on the production and formation space-time points extracted from the JETSET part of the PYTHIA. The ratio of the effective pre-hadron-nucleon cross section $\sigma_{\rm eff}$ to the usual total hadron-nucleon cross section $\sigma_{0}$ varies linearly with time between production and formation times as

[TABLE]

where $X_{0}=r_{\rm lead}a/Q^{2}$ is the pedestal value of the ratio, $r_{\rm lead}$ is the ratio of the number of leading quarks to the total number of quarks in a hadron, and $a=1$ GeV2 is a constant factor.

(ii) Quantum diffusion model (QDM) of ref. Farrar et al. (1988). This model has been applied for exclusive hard processes where, in the interaction point, the hadrons are in small-size configurations. In this case, the expression for the effective cross section is formally the same as Eq.(7). However, the production time $t_{\rm prod}$ is set to the time when the hard interaction happens, while the formation time is defined as $t_{\rm form}=t_{\rm prod}+l_{\rm form}/c$, where the formation length $l_{\rm form}$ is given by Eq.(6). In numerical calculations we set $a=0.56$ fm c/GeV in Eq.(6) as concluded from the analyses of pion electroproduction at JLab (see sec. 1). As far as the size of the quark-gluon configuration in the hard interaction point is concerned, in the case of DIS it is more difficult to match it to other observations. Since the model has a tendency to overestimate the neutron rate pretty significantly we will make a simplifying assumption that the interaction in the production point is strongly suppressed, that is $X_{0}=0$. As we will see below, even with this assumption, the model predicts a significantly higher neutron rate than observed experimentally.

(iii) A simple cutoff

[TABLE]

where $\Theta(x)$ is a Heaviside step function ($\Theta(x)=0$ for $x<0$ and $\Theta(x)=1$ for $x>0$). The cutoff momentum $p_{\rm cut}$ should be of the order of few GeV/c. It can be chosen from comparison with experimental data.

2.3 Particle propagation in the nucleus

The propagation of pre-hadrons and hadrons is described in the framework of the GiBUU model Buss et al. (2012). A pre-hadron is converted to an ordinary hadron at the time moment $t_{\rm form}$ which is determined for every produced pre-hadron. The mass and quantum numbers (spin, isospin, etc.) are not changed in that conversion. During formation time period, the pre-hadrons differ from respective hadrons only in the total interaction cross sections with the target nucleons and decay width (that disappears for the pre-hadrons). The branching ratios of different outgoing channels are the same for pre-hadronic and hadronic interactions. The nuclear mean field and Coulomb potentials are also supposed to be the same for pre-hadrons and hadrons 444It is of course questionable whether the pre-hadrons feel the same potentials. However, since the kinetic energy of a pre-hadron is much larger than the potential energy, the action of the mean field potential does not practically influence the pre-hadron propagation.. Bearing all that in mind we call ’hadron’ both hadronic and pre-hadronic states below in this subsection.

The GiBUU model explicitly solves the coupled set of kinetic equations for the system of hadrons in time and six-dimensional phase space. Between collisions and resonance decays the particles propagate according to the Hamiltonian equations of motion with the mean field potentials. The nuclear potential is described by the relativistic mean field model NL3 of ref. Lalazissis et al. (1997). The Coulomb potential acting on charge particles is also included. Particle collisions are simulated in the geometrical minimum distance scheme. The two particles approaching their minimal distance $d$ during a given time step will collide during this time step if $d<\sqrt{\sigma_{\rm tot}/\pi}$ where $\sigma_{\rm tot}$ is the total interaction cross section of these particles. Collisions are simulated by a Monte Carlo algorithm taking into account elastic and inelastic scattering channels. If the invariant energy $\sqrt{s}$ is larger than some threshold value (2.2 GeV for meson-baryon collision, 3.4 GeV for baryon-baryon collision, and 2.38 GeV for antibaryon-baryon collision), the hadron-hadron collision is simulated via the PYTHIA and (for antibaryon-baryon collisions only) the FRITIOF version 7.02 Andersson et al. (1993) models. The threshold values of $\sqrt{s}$ are selected to ensure the smooth behaviour of the hadronic cross sections over the transition from the resonance regime at low energies to the string-dominated regime at high energies (see sec. 3.3.2 of ref.Buss et al. (2012) and refs. Larionov et al. (2012); Weil (2013) for detail). The proton and neutron density profiles are chosen in a Woods-Saxon form with geometrical parameters taken from the Skyrme-Hartree-Fock systematics Lenske . The Fermi motion is taken into account in the local Fermi approximation. The Pauli blocking factors are included for the nucleons in the final state of scattering events. The nuclear mean field and Pauli blocking factors are supposed to be constant during the time evolution of the hadronic system. This allows to save CPU time for the elementary-particle-induced reactions on nuclei without significant accuracy lost.

There is also one important technical aspect which deserves to be mentioned. Most of the previous GiBUU calculations for the lepton-nucleus interactions were done with the so-called perturbative particles to increase statistics and save CPU time, see e.g. Appendix D.2 of ref. Buss et al. (2012). This method works very well for fast particles but not for slow ones, since the nucleus stays intact in calculations with perturbative particles. In particular, using the perturbative particles would lead to the drastic overestimation of the number of holes in the nucleus, see discussion below in sec. 2.4. Therefore, to avoid such problems, we apply in calculations the real-particles method which takes into account the trailing of the nucleus in the course of the cascade and thus provides much better description of the energy deposition in the nuclear residue. The importance of the trailing effect in proton-nucleus collisions at the energies of several GeV has already been found in early intranuclear cascade (INC) studies (cf. Barashenkov et al. (1974) and refs. therein).

2.4 Determination of an excited nuclear residue

The outgoing hadrons from a DIS event on the bound nucleon propagate through the nucleus inducing the cascade of interactions. This eventually leads to the emission of slow neutrons which can be produced either directly or by statistical decay of the excited residual nucleus.

The determination of the parameters of the residual nucleus can be done by various methods. In $pA$ collisions below pion production threshold, it is possible to apply the energy balance for the determination of the residue excitation energy Bertini (1963). However, this requires the utmost precision in the calculation of the total kinetic energy of emitted particles which renders this method impractical at the proton beam energy \raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}1 GeV. In the case of strongly violent nuclear dynamics, like central $pA$ collisions at GeV energies and central $AA$ collisions above Fermi energy, the nuclear excitation energy may reach several MeV/nucleon. The residual nuclear system strongly deviates from the ground state due to thermal excitation and may, for example, have a low density and/or a collective flow. In such situations, the only way is to use theoretical models with well defined ground state which allows to calculate directly the energy of the residual nuclear system and subtract its ground state energy by using the corresponding many-body Hamiltonian (cf. refs. Sangster et al. (1992); Botvina et al. (1995)) or relativistic energy density functional (cf. ref. Gaitanos et al. (2008)). However, if the system looses only relatively few ($\sim 10$) nucleons, the accuracy of such direct methods ($\sim 1$ MeV/nucleon) is not enough, while the old-fashioned residue determination in terms of particle-hole excitations works significantly better (cf. Eqs.(62),(63) in ref. Bondorf et al. (1995)).

Thus, we start adding up hole excitations in the nucleus neglecting, for simplicity, particle excitations. The corrections for the latter will be introduced at the last step. The mass number $A_{\rm res}$, charge number $Z_{\rm res}$, excitation energy $E^{*}_{\rm res}$, and momentum $\boldsymbol{p}_{\rm res}$ of the residual nucleus can be determined in every parallel ensemble555Parallel ensembles can be regarded as “events” in GiBUU simulations. by counting the hole excitations as

[TABLE]

where $A$($Z$) is the target nucleus mass (charge) number; $n_{\rm h}$ is the total number of holes; $Q_{i}=1(0)$ for the proton (neutron) hole; $E_{i}$ and $\boldsymbol{p}_{i}$ are, respectively, the single particle energy and momentum; $E_{F,i}$ is the Fermi energy. The dependence of the Fermi energy on the hole index “i” is mostly due to the difference of the Fermi energies for protons and neutrons666Since we use the empirical Woods-Saxon density profiles in combination with relativistic mean field potentials, the Fermi energies of protons and neutrons have also some spurious coordinate dependence. Comparison with the calculation using the relativistic Thomas-Fermi method for the nuclear ground states Gaitanos et al. (2010) shows, however, that the resulting excitation energy is practically insensitive to the details of the nuclear ground state density profile.. A hole is added to the nuclear residue for every two-body collision involving the target nucleon which did not collide yet, starting from the initial DIS event.

Let us discuss now the contribution of particle excitations which is neglected in Eqs.(9)-(12). The particle excitations are the bound nucleons with energies above Fermi level. A particle excitation is created if the recoil nucleon is not energetic enough and thus gets captured in the nuclear potential well. Hence, not every hole excitation leads to the reduction of the nucleus mass number. We have, therefore, redefined the quantities $A_{\rm res}$, $Z_{\rm res}$, and $\boldsymbol{p}_{\rm res}$ using the position space and counting all bound nucleons which are left in the target after the emission of the direct ones. The directly emitted particles have been determined by the requirement of their separation in the position space by at least 3 fm from all other nucleons. The choice of the minimum separation distance does not influence the result once this distance is larger than the internucleon distance of 1-2 fm and the time of the system evolution is large enough (see discussion in ref. Larionov et al. (2012)). If not specifically mentioned, the results reported below are obtained with the GiBUU calculation stopped at the time $t=100$ fm/c. The mass and charge numbers of the residue and its momentum are then replaced, respectively, by the total mass and charge numbers and momentum of the bound nucleons. The residue excitation energy of Eq.(11) is, however, not modified. This is because, in contrast to the INC models Barashenkov et al. (1972), we allow for collisions of the bound nucleons above Fermi level with other nucleons. As a result, the number of particle excitations grows, but they are approaching the Fermi level relaxing their kinetic energies to the low-lying particle-hole excitations.

2.5 Evaporation from the excited nuclear residue

The decay of excited nuclear residues has been described with a help of SMM Botvina et al. (1987); Bondorf et al. (1995). This model is focused on the description of the multiple simultaneous breakup of highly-excited nuclear systems (with excitation energy of the order of few MeV/nucleon). However, the SMM includes also the sequential decays of excited fragments which are treated similar to the Weisskopf model Weisskopf (1937).777In some contrast to the Weisskopf model, SMM includes also the decay width of an excited compound nucleus with respect to the emission of fragments heavier than $\alpha$-particle taking into account that the emitted fragment may be in an excited state stable to nucleon emission (see Eq. (34) in ref. Botvina et al. (1987)). To save CPU time, we have applied the SMM in the evaporation mode, switching-off the simultaneous multifragment breakup, since in the considered reactions the yield of nuclear residues with excitation energies exceeding $1-2$ MeV/nucleon is small (cf. Fig. 1 below). In some selected cases, we performed the full SMM calculation, including multifragmentation, and found out that this does not lead to noticeable changes of the neutron multiplicities and spectra.

3 Results

3.1 Muon-induced DIS on nuclei

In this section we analyze the E665 data on slow neutron production in inclusive $\mu^{-}+^{208}$Pb and $\mu^{-}+^{40}$Ca DIS processes at the beam energy of 470 GeV Adams et al. (1995a). The scattered muon was sampled as described in sec. 2.1 with cuts $\nu>20$ GeV, $Q^{2}>0.8$ GeV2. In the studied kinematics, a rather modest, on the scale of 20%, nuclear shadowing dominated by the resolved photon interactions with two nucleons is present for $R_{A}=\sigma(\mu Pb)/A\sigma(\mu N)$. In the setup of the E665 experiment, both inelastic and diffractive events were included. Applying Abramovsky-Gribov-Kancheli cutting rules Abramovsky et al. (1973) we observe that the events where projectile interacted with two nucleons and diffractive events contribute equal cross sections to the final state. While in the diffractive events no neutrons are produced, in the interactions where two nucleons are wounded roughly twice as many neutrons are produced. Overall, this results in the disappearance of shadowing for inclusive cross section (similar to the cancellation for central rapidities in ref. Abramovsky et al. (1973)) and in the enhancement of the average number of wounded nucleons per event:

[TABLE]

(see e.g. ref. Bertocchi and Treleani (1977)) increasing the average neutron multiplicity by roughly the same factor.

This implies that the constraints on the secondary interactions are somewhat stronger than those which we find in our analysis. For now we will neglect this effect.

Fig. 1 displays the mass number loss distribution (a), the excitation energy distribution (b) and the correlation between the excitation energy and the mass number loss (c) of the residual nucleus for the 208Pb target. Calculations in all models produce a sharp peak in the mass number loss distribution for the removal of only one nucleon and a broad tail towards larger mass number loss. For example, in the calculation with $p_{\rm cut}=1$ GeV/c, $57\%$ of residual nuclei have $A_{\rm res}=207$, while in the QDM calculation – $18\%$. In most of such events, a hole is produced due to initial hard DIS and the produced particles do not interact with the residue. Thus, the calculation with the strongest restriction on the momentum of interacting particle produces the largest fraction of one-hole events. The excitation energy distributions have a two-slope structure with a kink at $E^{*}_{\rm res}=50-70$ MeV. The region $E^{*}_{\rm res}<E_{F}\simeq 37$ MeV is saturated by one-hole excitations mostly. Higher excitation energies are only possible for the multiple (2,3,…) hole excitations. The dependence of the excitation energy on the number of removed nucleons is close to a linear one, $\langle E^{*}_{\rm res}\rangle\simeq 25~{}\mbox{MeV}(A-A_{\rm res})$.888Note that without contribution of particle excitations, this dependence would be $\langle E^{*}_{\rm res}\rangle\simeq 10~{}\mbox{MeV}(A-A_{\rm res})$. The difference arises due to the change of the number of removed nucleons, see discussion in sec. 2.4. Due to the Fermi motion, however, the distribution in the excitation energy for the fixed number of removed nucleons is quite broad. We checked that this distribution practically does not depend on the treatment of the pattern of hadron formation in calculations.

The spectrum of emitted neutrons is shown in Fig. 2a. Below 1 MeV, the spectrum is almost entirely due to statistical evaporation from the excited nuclear residue. In this low energy region, the number of neutrons depends weakly on the excitation energy of the nuclear residue (and thus on the treatment of hadron formation), since with decreasing $E^{*}_{\rm res}$ the energy spectrum of evaporated neutrons becomes more steep. The sensitivity of the spectrum to various treatments of hadron formation increases with neutron energy. The QDM and default calculations overestimate neutron yields. Only applying a very strong restriction on the momenta of interacting particles ($p<p_{\rm cut}=1$ GeV/c) allows us to describe the E665 data. Using the in-medium reduced low-energy $NN$ cross sections of refs. Li and Machleidt (1993, 1994) leads to $\sim 20\%$ reduction in the neutron multiplicity. This, however, would bring the calculation with $p_{\rm cut}=1$ GeV/c to even closer agreement with the E665 data for the lead target.

In Fig. 2b, the neutron spectra for different maximum times of the GiBUU calculation are shown. The dynamically emitted component of the spectra at low energies (E_{n}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}6 MeV) is moderately sensitive to variation of the maximum time, however, the full spectrum is practically stable.

We have also performed the calculation by setting $p_{\rm cut}$ to a very large value (1 TeV) which can be considered as no formation at all. This calculation produces almost indistinguishable with QDM yields of neutrons below 10 MeV (of course, at higher neutron energies the QDM calculation gives smaller yields than the calculation with instantly formed hadrons experiencing more final state interactions). However, the default GiBUU calculation gives smaller yields of slow neutrons as compared to the QDM one. It appears that the difference originates mainly from much larger expansion time, i.e. the string proper time between the hard interaction and first string breaking space-time points, in the GiBUU model Gallmeister and Falter (2005). In the nucleus rest frame, the expansion time corresponds to the time period $t<t_{\rm prod}$ when the effective cross section of Eq.(7) is set to zero in the default GiBUU calculation.

Having successfully described the E665 data for the lead target with $p_{\rm cut}=1$ GeV/c, we would expect that the description of the data for a lighter target nucleus would require even a lower value of the cutoff momentum. This is because a hadron likely passes through the nucleus without interaction, if its formation length is comparable with the nuclear size. Thus, we conclude from Eq.(6) that the cutoff momentum should be proportional to the nuclear radius. This naive picture is, however, not supported by the experimental data. Fig. 3 shows the average neutron multiplicity vs the photon energy $\nu$ for the calcium and lead targets. In the case of the 40Ca target, the neutron multiplicities are underestimated by a factor of 2.5 in the calculation with $p_{\rm cut}=1$ GeV/c and by 30% – in the QDM calculation. Thus, we fail to describe the E665 data for the calcium target based on the present-day formation length concept. This suggests (if the E665 neutron data are confirmed) a substantially more complicated mechanism of formation of slow hadrons probably involving a significant degree of coherence between the produced and spectator partons which ultimately form slow hadrons.

3.2 Quasireal photon interactions with nuclei at the LHC and RHIC

Nowadays, the feasibility of the studies of the heavy ion UPCs at colliders is widely realized. In such processes, one of the ions serves as a source of a quasireal photon which interacts with the second ion leaving the first ion intact. The UPCs are now intensively studied at the LHC. The most recent review of the UPC physics at the LHC can be found in ref. Baltz (2008). For a summary of the recent results and references we refer the reader to ref. Angerami (2018).

The maximum photon longitudinal momentum can be estimated in the c.m. frame of colliding nuclei (collider laboratory frame) from the condition that the photon wave length is equal to the radius of the Lorentz-contracted emitting nucleus Baltz (2008), i.e. $k_{L}^{\rm max}\simeq\gamma_{L}/R_{A}$, where $\gamma_{L}$ is the Lorentz factor. Assuming symmetric colliding system, we get $k^{\rm max}=\gamma_{L}2k_{L}^{\rm max}\simeq 2\gamma_{L}^{2}/R_{A}$ in the rest frame of another colliding nucleus. The estimated values of the maximum photon momentum $k^{\rm max}$ and of the $\gamma^{*}N$ c.m. energy $W$ are listed in Table 1 (see also Table 1 in ref. Baltz (2008))

In this paper we are discussing the simplest case when the photon experiences a hard interaction with one nucleon only and model it using the PYTHIA. In this model, the fragmentation of a nucleon weakly depends on the light cone momentum fraction $x_{N}$ of the nucleon carried by the removed parton for moderate $x_{N}$ and on the type of the removed parton. In principle, due to the color effects, the cases of the removal of a gluon (which is the dominating process in $AA$ UPCs) and a quark from the nucleus may differ. In particular, in the Lund model, the removal of a gluon, the $\gamma g\to\mbox{two jets}$ process, corresponds to the attachment of two strings to the nucleon leading to somewhat softer jets as compared to the quark case, $\gamma q\to q^{\prime}$. (This effect would be possible to test in $pA$ UPCs.) For now we neglect this difference and rely on the inclusive set of the PYTHIA events under certain kinematic conditions.

Since the PYTHIA model describes a $\gamma^{*}$-nucleon collision with the photon emitted by the scattered lepton, it is impossible to initialize an exactly real photon. However, by varying the kinematics of the lepton scattering one can change the photon virtuality $Q^{2}$ (or Bjorken $x=Q^{2}/2\nu m_{N}$) and the relative flux of the longitudinal photons $\epsilon$. We have checked that the variation of $\epsilon$ in the interval $0.02-0.10$ has practically no effect on the results. In the calculations below we set $\epsilon=0.05$. The $z$-axis is directed along the photon momentum.

In the dijet events measured at the LHC the jets are typically required to have $p_{t}>20$ GeV/c (cf. Aad et al. (2016)). For the back-to-back jets, this roughly corresponds to setting the minimum dijet invariant mass $M^{(0)}_{\rm dijet}=40$ GeV. We have chosen $W=100$ GeV and $W=500$ GeV as representative values of the $\gamma^{*}N$ c.m. energy at the LHC. The Bjorken $x$ has been set to $M^{(0)2}_{\rm dijet}/W^{2}$ which gives $x=0.16$ for $W=100$ GeV and $x=6.4\times 10^{-3}$ for $W=500$ GeV. This is equivalent to modeling the dijet production with invariant mass $M_{\rm dijet}>M^{(0)}_{\rm dijet}$ by the direct photon on a gluon carrying the light cone momentum fraction of the nucleon $x_{N}=M^{2}_{\rm dijet}/W^{2}$. The lower boundary imposed on the dijet invariant mass guaranties the smallness of the shadowing effects.

We have also performed calculations at RHIC kinematics with $W=30$ GeV, $x=0.16$. This corresponds to setting the threshold dijet invariant mass $M^{(0)}_{\rm dijet}=12$ GeV.

The mass number loss and excitation energy distributions of residual nuclei and the correlation between $A-A_{\rm res}$ and $E^{*}_{\rm res}$ in the LHC kinematics are almost identical to those for the E665 kinematics shown in Fig. 1. (Thus, the respective plots for the LHC are not shown.) The shape of the neutron multiplicity distributions shown in Fig. 4(a) for the LHC kinematics reflects the shape of the excitation energy distributions of the nuclear residues (Fig. 1b) with a broad maximum near $E^{*}_{\rm res}=0$ and a long tail of high energy excitations.

The variation of the photon kinematics and of the target nucleus has only a little effect on the neutron multiplicity distribution as demonstrated in Fig. 4(b) for $p_{\rm cut}=1$ GeV/c. This holds true also for other prescriptions for the hadron formation (not shown). Some excess of events with low neutron multiplicities is visible in the case of 197Au target with smaller $N/Z$ ratio.

In the ultraperipheral collision of heavy ions, it is convenient to study the transverse momentum spectrum of neutrons from the target nucleus since this spectrum is longitudinal-boost-invariant.

Experimentally, the neutrons are supposed to be detected in the zero degree calorimeter (ZDC) in the direction of the target nucleus. In the present exploratory work we do not intend to simulate all details of the experimental setups at the LHC and RHIC, which can be, in principle, done via the GiBUU+SMM framework.

The separation of neutrons in the direction of the target nucleus from those in the direction of the dijet can be done in terms of the longitudinal-boost-invariant Feynman variable $x_{F}$ defined as

[TABLE]

where $E$ ($E_{A}$) and $p^{z}$ ($p^{z}_{A}$) are the particle (target nucleus) energy and the longitudinal component of momentum, respectively. In the target nucleus rest frame the Feynman variable is expressed as

[TABLE]

The condition that the particle is emitted in the direction of the target nucleus in the collider laboratory frame is thus $x_{F}>x_{F}^{0}$, where $x_{F}^{0}=m_{t}/\sqrt{s_{NN}}$ with $m_{t}=\sqrt{m^{2}+p_{t}^{2}}$ being the particle transverse mass. In the LHC and RHIC kinematics the boundary values $x_{F}^{0}$ for the neutron with zero transverse momentum are $1.7\times 10^{-4}$ and $4.7\times 10^{-3}$, respectively. These values are still larger than those of $x_{F}^{\rm min}=m_{N}^{2}/W^{2}$ as obtained from the extreme kinematics of the $\gamma^{*}N\to N$ transition.

In concrete calculations we choose the cut $x_{F}>0.1$ relevant, for example, for the ZDC of ATLAS detector White (2011).

As we see from Fig. 5, the chosen $x_{F}$-cut practically does not influence the low-$p_{t}$ neutron spectra originating mostly from target fragmentation.

At the intermediate transverse momenta, $p_{t}\sim 1$ GeV/c, the $x_{F}$-cut reduces the neutron yield somewhat. The reduction effect, however, depends on the scenario of hadron formation: The neutrons with transverse momenta of a few GeV/c originate mostly from the nucleon fragmentation in the reaction $\gamma+N\to\mbox{two jets}+n+X$ on the bound nucleon. In the case of QDM calculation, such neutrons may still interact with surrounding nucleons and, therefore, increase their $x_{F}$ values due to deceleration. In the case of calculation with $p_{\rm cut}=1$ GeV/c, however, these neutrons do not interact with nucleus, and thus a part of them remains in the ’invisible’ region $x_{F}<0.1$.

The measurement of the neutron production at p_{t}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}1 GeV/c offers an additional test of the hadron formation scenario. Namely, if only nucleons with momenta below 1 GeV/c interact with the residual nucleus, the inclusive cross section in this kinematics should be linear in $A$ reflecting the direct mechanism of the neutron production 999Observing this effect would require a highly segmented ZDC. The ZDCs used in the current collider experiments detect with a very high efficiency neutrons produced with momenta below the Fermi momentum in the nucleus rest frame. To detect neutrons originating from the $\gamma N$ reaction on the bound nucleon one would need a different detector like the one which is employed by the LHCf Adriani et al. (2018)..

Fig. 6(a) shows the neutron low-$p_{t}$ spectra for different prescriptions for hadron formation. At very small transverse momenta, statistical evaporation hides the effect of hadron formation. However, already for $p_{t}$ above $\sim 100$ MeV/c we observe a much stronger dependence of the neutron emission on the pattern of hadron formation. As shown in Fig. 6(c), the spectra at low $p_{t}$’s are almost independent on $W$, $x$ and on the choice of the nuclear target.

This means that the convolution with the photon flux should not influence the low-$p_{t}$ spectra of neutrons.

The above behaviour of the neutron $p_{t}$ spectra continues also towards higher transverse momenta as shown in Figs. 6b,d, if the cutting condition $x_{F}>0.1$ is applied. However, the $p_{t}$ region above 1 GeV/c is significantly influenced by the $x_{F}$ cut. In particular, as we can see from Fig. 5, in full spectra, the dependence on the hadron formation scenario tends to disappear with increasing $p_{t}$. On the other hand, the full neutron spectra at p_{t}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}1 GeV/c become influenced by photon kinematics (not shown). The detailed discussion of this transverse momentum range goes beyond the scope of our work. We would only mention one feature which exists both in the full and $x_{F}$-cut spectra at $p_{t}$ above 1 GeV/c. Namely, the per-event spectra on the proton and nucleus are similar at high $p_{t}$’s, however, the latter spectra are significantly enhanced at low $p_{t}$’s (cf. Fig. 7a and Fig. 6d). This is a clear indication of the production mechanism of low-$p_{t}$ neutrons in secondary interactions. Another indication is the sensitivity of the low-transverse-momentum neutron yield to the hadron formation.

3.3 Discussion

To understand better the sources of slow neutrons we present in Fig. 7b the decomposition of the neutron $p_{t}$-spectrum in $\gamma^{*}+^{208}$Pb deep inelastic collisions in partial components. The direct production channel (denoted ’w/o FSI’) includes neutrons produced in a primary $\gamma^{*}N$ collision on a bound nucleon and emitted without final state interactions. The partial components of the secondary interactions are defined in the same way as in $pA$ collisions (see Fig. 10c). The high-$p_{t}$ part is dominated by the direct neutrons and by the neutrons produced in $\Delta$ decays, while the low-$p_{t}$ part – by the neutrons scattered elastically on a nucleon at the end of their interaction chain. Similar to the case of $pA$ collisions, the neutrons produced in inelastic channels and resonance decays and then scattered elastically constitute $\sim 15\%$ of the latter component. Therefore, most of the low-$p_{t}$ neutrons are produced in the multiple $NN$ elastic scattering processes.

Let us finally take a more detailed look at the primordial particles produced in the elementary $\gamma^{*}+p$ deep inelastic collision. Fig. 8 shows the inclusive $x_{F}$- and momentum-distributions of pions, protons, and neutrons. The condition $x_{F}<1-x$ is fulfilled, which follows from the kinematics of the target fragmentation. Overall, this leads to the enhancement of the slow particle production with decreasing $x$, which is most visible for protons (Fig. 8(c),(d)).

The directly produced slow neutrons are even more sensitive to the photon kinematics (Fig. 8(e),(f)). Like in the proton case, the neutron spectrum also shifts to smaller momenta with decreasing $x$, but in addition gets significantly suppressed because the $p\to n$ transition involves the valence quark. We can see from Fig. 7b that primordial protons and neutrons with momenta \raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}1 GeV/c govern the production of low-$p_{t}$ neutrons. While the momentum spectra of primordial protons and neutrons are influenced by the photon kinematics, it is not the case for production of the low-$p_{t}$ neutrons. This nontrivial behaviour can be explained by the fact that the sum of the proton and neutron primordial spectra above Fermi momentum and below 1 GeV/c is weakly influenced by the photon kinematics.

We see that pions are produced most abundantly (Fig. 8(a),(b)). The yields of protons and neutrons are suppressed by an order of magnitude relative to the pion yield. In the case of nuclear target, pions with momentum $\sim 0.3$ GeV/c are strongly absorbed in the processes $\pi N\to\Delta$. This explains large $\Delta$-decay and $N\Delta\to NN$ absorption components in the neutron transverse momentum spectrum below 0.5 GeV/c in Fig. 7b.

4 Directions for further studies

We have considered the elementary process most close to the leading twist DIS – the production of dijets in the direct photon – gluon interaction in the kinematics where nuclear shadowing effects are small. Here we outline several other interesting processes.

(a) In the model we ignored possible differences between the final state interactions for the case of the gluon and quark removals. One possible difference arises in the Lund string model where effectively two strings are approximating an octet residual system. Other effect is a larger energy loss in the propagation of a color octet dijet through the media leading to a higher nuclear excitation energies due to larger number of wounded nucleons. Although the production of slow neutrons is expected to be dominated by the interactions of products of the nucleon fragmentation with surrounding nucleons, the jet propagation may still provide significant contribution since there should be no rapidity gap between the products of the dijet- and nucleus fragmentation Angerami (2018).

(b) New effects are expected in the kinematics where a dijet is produced in the resolved photon – nucleus interaction. In this case, the active parton belonging to the photon wave function carries the light cone fraction $x_{\gamma}$ of the photon momentum which is significantly smaller than one. The selection of a parton with a smaller $x_{\gamma}$ in the photon corresponds to the selection of configurations with a larger number of partons in the photon wave function. For x_{\gamma}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}0.1, these configurations should interact with the target with the strength comparable to the strength of interaction of configurations responsible for nuclear shadowing in coherent diffractive photoproduction of $\rho$-mesons off nuclei Frankfurt et al. (2016). Based on this analysis we expect that the average number of wounded nucleons should be $n_{\rm wound}\sim 3$. Correspondingly, the average number of emitted neutrons should increase by a factor of three, since, to a first approximation, wounded nucleons fragment independently.

(c) One can consider minimum bias $\gamma A$ collisions at $W\sim 100-500$ GeV. Based on the analysis of the color fluctuations in photon in ref. Alvioli et al. (2017) we expect that in these collisions the average number of wounded nucleons would be large and fluctuating strongly. The number of emitted neutrons should be correlated with $p_{t}$ of the leading particles produced in the photon fragmentation region (smaller $p_{t}$ – larger neutron number).

(d) It appears feasible to detect charm particles in the UPCs. At large $W\geq 100$ GeV and small enough p_{t}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}\mbox{few GeV/c}, one reaches the kinematics where the nuclear shadowing becomes significant (a factor of two reduction of the cross section as compared to the impulse approximation) resulting in wounding in average two nucleons and hence increasing the number of neutrons by a factor of two.

5 Summary and conclusions

We have performed hybrid, i.e. dynamical plus statistical, calculations of slow neutron production in hard virtual photon induced inelastic collisions with nuclei. For the dynamical stage, the quantum-kinetic GiBUU transport model Buss et al. (2012) has been used. The statistical stage has been described by using the SMM model Botvina et al. (1987); Bondorf et al. (1995). For the calculation of the characteristics ($A_{\rm res},Z_{\rm res},E^{*}_{\rm res},\boldsymbol{p}_{\rm res}$) of the nuclear residue, a procedure based on the particle-hole excitations in the Fermi sea has been applied. We have tested these hybrid calculations by comparison with available data on the neutron production in muon- and proton-nucleus interactions.

The multiplicity of slow ($E<10$ MeV) neutrons in $\mu^{-}$ DIS on Pb measured by the E665 collaboration Adams et al. (1995a) is found to be a factor of two smaller than what is expected from GiBUU+SMM calculations both with the default treatment of hadron formation (based on hadron production and formation space-time points from JETSET) and with the quantum diffusion model of the color transparency. The neutron multiplicity and the energy spectrum in $\mu^{-}+^{208}$Pb DIS can be reproduced only if produced hadrons with momenta above 1 GeV/c are not allowed to interact with the rest of the nucleus.

This indicates the presence of a novel dynamics in the production of hadrons in the nuclear fragmentation region which cannot be reproduced by adjusting the parameters of standard models. In hindsight this may be not too surprising as we are dealing with particles produced in the target fragmentation region, while the models were designed to describe the space-time picture of hadron formation in the current fragmentation region. Hence, further studies are necessary to get a better understanding of the hadron formation in the target fragmentation region.

In the long run, it will be possible to perform such studies at the EIC by measuring the $A$-dependence of hadron spectra in the nucleus fragmentation region. However, already now the study of hard $\gamma A$ collisions in heavy ion UPCs at the LHC and RHIC with detection of neutrons in the ZDC would allow to clarify the space-time picture of hadron formation in the nucleus fragmentation region. To this end, we have performed exploratory studies of $\gamma$ + Pb(Au) interactions in the direct photon regime in the LHC (RHIC) kinematics. We have found a strong dependence of the neutron multiplicities and $p_{t}$-spectra on the model of hadron formation indicating that the slow neutron production could be a sensitive tool for studying the space-time picture of formation of slow hadrons in hard processes. We have also listed several directions for further studies of the interplay of hard and soft dynamics in the UPCs. In particular, we predict a strong increase of the neutron multiplicity in the resolved photon interactions as compared to the direct photon interactions.

All modifications of the GiBUU code done in the course of this work concerning the hadron formation time and determination of the nuclear residue are included in the public release GiBUU 2019. Our hybrid GiBUU+SMM calculations can be readily adopted for the geometry of the ZDC in particular experiments.

Acknowledgements.

The support by the Frankfurt Center for Scientific Computing is gratefully acknowledged. We thank Marcus Bleicher, Alexander Botvina, Leonid Frankfurt, Kai Gallmeister, and Ulrich Mosel for illuminating discussions. We thank A. Angerami for discussion of the experimental feasibility of the studies proposed in this paper and Yu. Dokshitzer for discussion of various mechanisms of the final state interactions. The authors are especially grateful to Ulrich Mosel who read the manuscript before publication and suggested substantial modifications of the abstract, introduction and conclusions. AL is grateful to the Department of Physics at the Pennsylvania State University, where this work was started, for the hospitality and financial support. AL also acknowledges partial financial support of HIC for FAIR within the framework of the Hessian LOEWE program. The research of MS was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award No. DE-FG02-93ER40771.

Appendix A $pA$ collisions

The purpose of this Appendix is to test our model set-up for the hadron propagation and statistical decay of the residual nuclei by comparison with experimental data on neutron production in proton-nucleus interactions at the beam energy about 1 GeV. In this intermediate energy range, the effects of hadron formation are unimportant. The proton-nucleus reaction proceeds as a cascade of elastic and inelastic (mainly, $\Delta$-resonance production and decay) processes and pion production and reabsorption. The late stage of the reaction, when the fast cascade particles have left the nucleus and the slow deexcitation of the nuclear residue starts – should be, however, similar to the reactions at high energies. We will consider the cases of the heavy targets, i.e. gold and lead. The proton was initialized at the distance of 0.5 fm from nuclear surface at the impact parameters $b=0.5,1.0,\ldots,7$ fm with 2000 events for every impact parameter. The calculated observables are weighted with impact parameter.

It is instructive, first, to test our model setup by comparison with previous calculations. We selected the INC calculations of ref. Barashenkov et al. (1974).101010See line “b” in Figs. 12 and 13 of ref. Barashenkov et al. (1974), corresponding to events without fission. Note that the fission cross section is \raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}10\% of the inelastic cross section for $p+$Au collisions at $E_{\rm lab}=1$ GeV. Fig. 9 shows the mass number loss– and excitation energy distributions of the nuclear residue for $p+^{197}$Au collisions at the beam energy of 1 GeV. It is evident that the correction of the residue mass number due to the particle excitations is quite large. Our corrected results are in a good agreement with INC results. Thus, we have included the correction for the particle excitations in all calculations discussed in sec. 3 and below in this Appendix. However, as we checked, this correction has practically no effect, since the statistical evaporation of neutrons is mostly governed by the excitation energy of the nucleus and not by its mass and charge numbers.

Fig. 10a,b show the neutron energy spectra at the laboratory angle of $120\degree$ for the proton-lead collisions at 1.4 and 2 GeV/c. The GiBUU calculation was run until 100 and 200 fm/c. Then the target residues were determined and used as an input to the SMM calculation which was run 10 times for every GiBUU event. Evaporated neutrons were included in the total spectrum of emitted neutrons. As we see the evaporation substantially increases the yield of low-energy (E\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}30 MeV) neutrons and explains the evident two-slope structure of the experimental neutron spectra. This is not surprising since the same data were explained earlier with similar approaches in ref. Strikman et al. (1999) and recently in ref. Galoyan et al. (2015).

Fig. 10c shows the neutron energy spectrum from 1.4 GeV/c $p+^{208}$Pb collisions with partial contributions of the $NN$ collisions, $N\Delta$ collisions, $NN\pi$ (pion absorption by nucleon pair) and $N\pi$ collisions (non-resonant channels), and $\Delta$ resonance decays. The partial contributions are selected according to the last process in which the neutron was finally produced. It means that, for example, a neutron may appear in the final state of an inelastic scattering or a resonance decay and afterwards scatter elastically. Then this neutron belongs to the $NN$ partial contribution. We estimated that such neutrons constitute $\sim 15\%$ of the $NN$ component. Thus, the main contribution to the spectrum of preequilibrium neutrons is given by $NN$ elastic collisions. The $N\Delta$ and $NN\pi$ collisions make comparable in the magnitude contribution above 100 MeV neutron energy. Other neutron production channels are less important.

By performing an additional calculation keeping only elastic scatterings we have established that the main contribution to the neutron spectrum comes from multiple scattering (from two- to five-step) processes as displayed in Fig. 10d. This suggests that the residual nuclear system is approaching thermally equilibrated state while dynamically emitting neutrons. This is a type of the situation where it is natural to apply statistical approaches.