1 Introduction

Measurements of particle production at very small polar angles with respect to the proton beam direction (forward direction) in positron-proton collisions are important inputs for the theoretical understanding of proton fragmentation mechanisms. Forward particle measurements are also valuable for high energy cosmic ray experiments, as they provide important new constraints for high energy air shower models [1, 2].

The H1 and ZEUS experiments at the \(ep\) collider HERA have studied the production of forward baryons (protons and neutrons) and photons, which carry a large fraction of the longitudinal momentum of the incoming proton [38]. These analyses have demonstrated that models of deep-inelastic scattering (DIS) are able to reproduce the forward baryon measurements if contributions from different production mechanisms are considered, such as string fragmentation, pion exchange, diffractive dissociation and elastic scattering of the proton [6, 7]. The forward photon production rate, however, is overestimated by the models by \(50\) to \(70~\)% [8]. The measurements also confirm the hypothesis of limiting fragmentation [9, 10], according to which, in the high-energy limit, the cross section for the inclusive production of particles in the target fragmentation region is independent of the incident projectile energy.

Measurements in the DIS regime provide a possibility to investigate the process at different centre-of-mass (CM) energies of the virtual photon-proton system, \(W\), within the same experiment. The studies of the energy dependence of particle production allow a test of the Feynman scaling [11] hypothesis, according to which particle production is expected to show a scaling behaviour, i.e. independence of the CM energy in terms of the Feynman-\(x\) variable, \(x_F=2p_{||}^*/W\). Here \(p_{||}^*\) is the longitudinal momentum of the particle in the virtual photon–proton CM frame with respect to the direction of the beam proton. In several previous measurements Feynman scaling was found to be violated in the fragmentation process in the central rapidity region [1221]. On the contrary, no sizable violation of Feynman scaling has been observed in the target fragmentation region in \(pp\) and \(p\bar{p}\) collisions by comparing the \(\pi ^0\) production cross sections at the SPS collider [22] with \(\pi ^{\pm }\) measurements at the ISR [2326]. However, these conclusions are debated [27] and the scarcity of other experimental forward particle production data motivates further studies of forward particle production.

In this paper the production of forward neutrons and photons in DIS is studied as a function of \(x_F\) Footnote 1 and \(W\). This is the first direct experimental test of Feynman scaling for photons and neutrons produced in the very forward direction. Predictions from different DIS and different cosmic ray (CR) hadronic interaction Monte Carlo (MC) models are compared to the results. The simultaneous measurement of forward neutrons and photons provides a useful input for MC model development also because of the respective different production mechanisms: forward photons almost exclusively originate from decays of neutral mesons produced in the fragmentation of the proton remnant (Fig. 1a), while forward neutrons are produced also via a colour singlet exchange process (Fig. 1b).

Fig. 1
figure 1

a Generic diagram for forward photon or neutron production \(ep\rightarrow e'\gamma X\), \(ep\rightarrow e'nX\) in deep-inelastic scattering. b Diagram of forward neutron production via pion exchange

The neutrons and photons studied here are produced at polar angles below \(0.75\) mrad and are measured in the Forward Neutron Calorimeter (FNC) of the H1 detector. The data used in this analysis were collected with the H1 detector at HERA in the years 2006 and 2007 and correspond to an integrated luminosity of \(131~ \mathrm pb^{-1}\). During this period HERA collided positrons and protons with energies of \(E_e=27.6~\mathrm {GeV}\) and \(E_p=920~\mathrm {GeV}\), respectively, corresponding to a centre-of-mass energy of \(\sqrt{s}=319~\mathrm {GeV}\).

2 Experimental procedure and data analysis

2.1 H1 main detector

A detailed description of the H1 detector can be found elsewhere [2833]. Only the detector components relevant to this analysis are briefly described here. The origin of the right-handed H1 coordinate system is the nominal \(e^+p\) interaction point. The direction of the proton beam defines the positive \(z\)-axis; the polar angle \(\theta \) is measured with respect to this axis. Transverse momenta are measured in the \(x\)\(y\) plane. The pseudorapidity is defined by \(\eta = -\ln {[\tan (\theta /2)]}\) and is measured in the laboratory frame.

The interaction region is surrounded by a two-layer silicon strip detector and large concentric drift chambers, operated inside a \(1.16\) T solenoidal magnetic field. Charged particle momenta are measured in the angular range \(15^\circ <\theta <165^\circ \). The forward tracking detector is used to supplement track reconstruction in the region \(7^\circ <\theta <30^\circ \) and improves the hadronic final state reconstruction of forward going low momentum particles. The tracking system is surrounded by a finely segmented liquid argon (LAr) calorimeter, which covers the polar angle range \(4^\circ <\theta <154^\circ \) with full azimuthal acceptance. The LAr calorimeter consists of an electromagnetic section with lead absorber and a hadronic section with steel absorber. The total depth of the LAr calorimeter ranges from \(4.5\) to \(8\) hadronic interaction lengths \(\lambda \). The absolute electromagnetic energy scale is known with a precision of \(1~\%\), while the absolute hadronic energy scale is known for the present data with a precision of \(4~\%\).

The backward region (\(153^\circ <\theta <177.8^\circ \)) is covered by a lead/scintillating-fibre calorimeter called the SpaCal; its main purpose is the detection of scattered positrons. The polar angle of the positron is measured with a precision of \(1\) mrad. The energy resolution for positrons is \(\sigma (E)/E\approx 7.1~\%/\sqrt{E[\mathrm {GeV}]}\oplus 1~\%\) [34] and the energy scale uncertainty is less than \(1~\%\). The hadronic energy scale in the SpaCal is known with a precision of \(7~\%\).

The luminosity is determined from the rate of the elastic QED Compton process with the electron and the photon detected in the SpaCal calorimeter, and the rate of DIS events measured in the SpaCal calorimeter [35].

2.2 Forward detector for neutral particles

Neutral particles produced at very small polar angles with respect to the proton beam direction can be detected in the FNC, which is situated at a polar angle of \(0^\circ \) at \(106\) m from the interaction point. A detailed description of the detector is given in [7, 8]. The FNC is a lead–scintillator sandwich calorimeter. It consists of two longitudinal sections: the Preshower Calorimeter with a length corresponding to about \(60\) radiation lengths or \(1.6 \lambda \) and the Main Calorimeter with a total length of \(8.9 \lambda \). The acceptance of the FNC is defined by the aperture of the HERA beam-line magnets and is limited to scattering angles of \(\theta \lesssim 0.8\) mrad with approximately \(30~\%\) azimuthal coverage.

The longitudinal segmentation of the FNC allows an efficient discrimination of photons from neutrons. Photons are absorbed completely in the Preshower Calorimeter, while neutrons have a significant fraction of their energy deposited in the Main Calorimeter. Therefore, energy deposits in the FNC, which are contained in the Preshower Calorimeter with no energy deposits above the noise level in the Main Calorimeter, are classified as electromagnetic clusters. According to the Monte Carlo simulation about \(98~\%\) of all reconstructed photon and neutron candidates originate from generated photons and neutrons, respectively. Due to the relatively large size of the FNC readout modules in combination with the small geometrical acceptance window, two or more particles entering the FNC are reconstructed as a single cluster. In the MC simulation about \(1~\%\) of all hadronic clusters in the FNC associated with neutrons are overlapped with a photon, which was scattered within the FNC acceptance together with the neutron. At lower energies the electromagnetic clusters reconstructed in the FNC mainly originate from single photons. At higher measured energies there is a significant contribution from two photons, with the fraction of two-photon events increasing from \(0.5~\%\) at \(100\) GeV to \(10~\%\) at about \(450\) GeV and to \(80~\%\) at \(900\) GeV. The two photons typically originate from the decay of the same meson.

The absolute electromagnetic and hadronic energy scales of the FNC are known to \(5~\%\) [8] and \(2~\%\) [7] precision, respectively. The energy resolution of the FNC calorimeter for electromagnetic showers is \(\sigma (E)/E \approx 20~\%/\sqrt{E~[\mathrm GeV]}\ \oplus 2~\%\) and for hadronic showers \(\sigma (E)/E \approx 63~\%/\sqrt{E~[\mathrm GeV]} \oplus 3~\%\), as determined in test beam measurements [36]. The spatial resolution is \(\sigma (x,y)\approx 10~\mathrm cm/\sqrt{E~[\mathrm GeV]} \oplus 0.6~\mathrm cm\) for hadronic showers starting in the Main Calorimeter. A better spatial resolution of about \(2\mathrm ~ mm\) is achieved for the electromagnetic showers and for those hadronic showers which start in the Preshower Calorimeter.

2.3 Cross section definition

The kinematics of semi-inclusive forward photon and neutron production are shown in Fig. 1a, where the four-vectors of the incoming and outgoing particles and of the exchanged virtual photon \(\gamma ^*\) are indicated. This measurement is restricted to the DIS kinematic range, determined by the photon virtuality \(6<Q^2<100~\mathrm {GeV}^2\) and inelasticity \(0.05<y<0.6\). They are defined as

$$\begin{aligned} Q^2=-q^2, \quad y=\frac{p \cdot q}{p \cdot k}\;, \end{aligned}$$
(1)

where \(p\), \(k\) and \(q\) are the four-momenta of the incident proton, the incident positron and the virtual photon, respectively. The CM energy of the virtual photon-proton system, \(W\), is related to \(Q^2\) and \(y\) as \(W \approx \sqrt{ys-Q^2}\), where \(s\) is the squared total CM energy of the positron-proton system. The present analysis is restricted to the range \(70<W<245~\mathrm GeV\).

The analysis is performed in the pseudorapidity range \(\eta >7.9\) for forward neutrons and photons. The pseudorapidity range \(\eta >7.9\) corresponds to polar angles \(\theta <0.75\) mrad, the geometrical acceptance limit of the FNC. In the virtual photon-proton CM frame the \(x_F\) of the neutron and photon are restricted to \(0.1<x_F<0.94\) and \(0.1<x_F<0.7\), respectively. These requirements limit the transverse momentum of neutrons to \(0<p_T^*<0.6~\mathrm GeV\) and of photons to \(0<p_T^*<0.4~\mathrm GeV\). The requirement that \(x_F\) is below \(0.7\) for photons ensures that the electromagnetic clusters reconstructed in the FNC mainly originate from single photons, according to MC predictions.

The kinematic phase space of the measurements is summarised in Table 1. Cross sections of neutrons and photons produced in the forward direction, normalised to the inclusive DIS cross section, \( 1/\sigma _{DIS}~\mathrm{d}\sigma /\mathrm{d}x_F\), are determined differentially in \(x_F\) in three ranges of \(W\). In addition, the cross section ratios integrated over \(x_F\), \( \sigma ^{\gamma ,n}_{DIS}/\sigma _{ DIS}\), are measured as a function of \(W\).

Table 1 Definition of the kinematic phase space of the measurements

2.4 Event selection

The data selection and analysis procedures are similar to those described in previous publications using the FNC [7, 8]. The data sample of this analysis was collected using triggers which require the scattered positron to be measured in the SpaCal. The trigger efficiency is about 96 % for the analysis phase space as determined from data using an independently triggered data sample. The selection of DIS events is based on the identification of the scattered positron as the most energetic, isolated compact calorimetric deposit in the SpaCal with an energy \(E_e'>11~\mathrm {GeV}\) and a polar angle \(156^\circ <\theta _e^{\prime }<175^\circ \). The \(z\)-coordinate of the primary event vertex is required to be within \(\pm 35\) cm of the nominal interaction point. The hadronic final state is reconstructed using an energy flow algorithm which combines charged particles measured in the trackers with information from the SpaCal and LAr calorimeters [37, 38]. To suppress events with hard initial state QED radiation, as well as events originating from non-\(ep\) interactions, the quantity \(\sum {E - p_z}\), summed over all reconstructed particles including the positron, is required to lie between 35 and \(70~\mathrm {GeV}\). This cut also efficiently removes events from photoproduction processes, where the positron is scattered into the backward beam-pipe and a particle from the hadronic final state fakes the positron signature in the SpaCal. The kinematic variables \(Q^2\) and \(y\) are reconstructed using a technique which optimises the resolution throughout the measured \(y\) range, exploiting information from both the scattered positron and the hadronic final state [39]. Events are restricted to the range \(6<Q^2<100~\mathrm {GeV}^2\) and \(0.05<y<0.6\). The DIS data sample contains about 9.3 million events.

A subsample of events containing forward photons or neutrons is selected by requiring either an electromagnetic or a hadronic cluster in the FNC with a pseudorapidity above 7.9 and an energy above \(92~\mathrm GeV\). The data sample, called ‘FNC sample’ in the following, contains about 83,000 events with photons and 230,000 events with neutrons.

2.5 Monte Carlo simulations and corrections to the data

Monte Carlo simulations are used to correct the data for the effects of detector acceptance, inefficiencies, QED radiation from the positron and migrations between measurement bins due to the finite detector resolution. All generated events are passed through a GEANT3 [40] based simulation of the H1 apparatus and are subject to the same reconstruction and analysis chain as the data.

The DJANGOH [41] program is used to generate inclusive DIS events. It is based on leading order electroweak cross sections and takes into account QCD effects up to order \(\alpha _s\). Higher order QCD effects are simulated using leading log parton showers as implemented in LEPTO  [42], or using the Colour Dipole Model (CDM) as implemented in ARIADNE [43]. Subsequent hadronisation effects are modelled using the Lund string fragmentation model as implemented in JETSET [44, 45]. Higher order electroweak processes are simulated using an interface to HERACLES [46]. The LEPTO program optionally includes the simulation of soft colour interactions (SCI) [47], in which the production of diffraction-like configurations is enhanced via non-perturbative colour rearrangements between the outgoing partons. The SCI option in LEPTO is used for the simulation of forward photons but not for neutrons because it has been shown [7] that the description of the forward neutron data is poor in this case. For the DJANGOH MC simulations the H1PDF 2009 parameterisation [48] of the parton distributions in the proton is used. In the following, the DJANGOH predictions based on LEPTO and ARIADNE are denoted as LEPTO and CDM, respectively. In all DJANGOH simulations forward particles originate exclusively from the hadronisation of the proton remnant and forward photons are therefore mainly produced from the decay of \(\pi ^0\) mesons.

RAPGAP [49] is a general purpose event generator for inclusive and diffractive \(ep\) interactions. Higher order QCD effects are simulated using parton showers and the final state hadrons are obtained via the Lund string model. As in DJANGOH higher order electroweak processes are simulated using an interface to HERACLES [46]. In the version denoted below as RAPGAP-\(\pi \), the program simulates exclusively the scattering of virtual or real photons off an exchanged pion (Fig. 1b). In this model the cross section for \(ep\) scattering to the final state \(nX\) takes the form

$$\begin{aligned} {\mathrm d}\sigma (e p\rightarrow e'nX) =f_{\pi ^+/p}(x_L,t)\cdot {\mathrm d}\sigma (e\pi ^+\rightarrow e'X). \end{aligned}$$
(2)

Here \(x_L\) is the longitudinal momentum fraction and \(t\) is the squared four-momentum transfer between the incident proton and the final state neutron; \(f_{\pi ^+/p}(x_L,t)\) represents the pion flux associated with the splitting of a proton into a \(\pi ^+ n\) system and \({\mathrm d}\sigma (e\pi ^+\rightarrow e'X)\) is the cross section of the positron-pion interaction. There are several parameterisations of the pion flux [5054] and the one used here is taken from [51]. The details of the pion flux parameterisation are described in [7]. Using other parameterisations of the pion flux affects mainly the absolute normalisation by up to \(30~\%\).

As was shown in [7], the best description of the forward neutron data is achieved by a combination of events with neutrons originating from pion exchange, as simulated by RAPGAP-\(\pi \), and events with neutrons from proton remnant fragmentation, simulated by DJANGOH without SCI. RAPGAP-\(\pi \) mainly contributes at high neutron energies, while DJANGOH is significant at low energies. In [7] the contributions of RAPGAP-\(\pi \) and CDM were added using weighting factors \(0.65\) and \(1.2\) for the respective predictions. In the present analysis the best description of the neutron energy distribution is obtained by the combination of RAPGAP-\(\pi \) and CDM using the respective weights \(0.6\) and \(1.4\), or by the combination of RAPGAP-\(\pi \) and LEPTO using the respective weights \(0.6\) and \(0.7\). The difference between the weighting factors for the combination of RAPGAP-\(\pi \) and CDM in this analysis and in [7] is due to the different neutron energy range and the resulting different neutron energy dependence in the two analyses. Compared to [7] the current analysis is extended to much lower neutron energies, at which the contribution from the fragmentation model dominates.

The measurements are also compared to predictions of several hadronic interaction models which are commonly used for the simulation of cosmic ray air shower cascades: EPOS LHC [55, 56], QGSJET 01 [57, 58], QGSJET II-04  [59, 60] and SIBYLL 2.1 [61, 62]. These models are based on Regge theory [63], on the Reggeon calculus of Gribov [64] and on perturbative QCD. They use an unitarisation procedure to reconstruct amplitudes for exclusive processes and to determine the total and elastic cross sections. Central elements of these models are the production of mini-jets and the formation of colour strings that fragment into hadrons. Whereas the Regge-Gribov approximation is applied to hadrons as interacting objects in the case of QGSJET and SIBYLL, it is extended to include partonic constituents in EPOS LHC. Compared to the earlier EPOS simulation [55], which was used for comparison with the previous H1 forward photon analysis [8], the new EPOS LHC model [56] includes a modified treatment of central diffraction and the diffractive remnant in order to reproduce rapidity gap measurements at the LHC. The CR models also differ in the treatment of saturation effects at high parton densities at small Bjorken-\(x\) and in the treatment of the hadronic remnants in collisions. The programs are interfaced with the PHOJET program [65] for the generation of the \(ep\) scattering kinematics. It was pointed out [66] that the hadronic interaction models have been developed for hadron-hadron interactions and therefore the simulation of DIS events might be affected by the superfluous contribution of multi-parton interactions. In order to investigate this assumption, the QGSJET 01 model has been modified [67] to exclude the multi-parton interactions. In the comparison with the measurements this modified model is denoted as ‘QGSJET 01 (no mi)’.

The measured distributions may contain background arising from several sources. The background from photoproduction processes is estimated using the PHOJET MC generator. It is found to be about \(1~\%\) on average and is subtracted from the data distributions bin-by-bin. Background from misidentification of photons or neutrons in the FNC is estimated from the DJANGOH  MC simulation to be \(2~\%\) on average and is subtracted from the data distributions bin-by-bin. Background also arises from a random coincidence of DIS events, causing activity in the central detector, with a beam-related background signal in the FNC, produced from the interaction of another beam proton with a positron or with residual gas in the beampipe. This contribution is estimated by combining DIS events with FNC clusters originating from interactions in the bunch-crossings adjacent to the bunch-crossings of the DIS events. It is found to be smaller than \(1~\%\) and is neglected.

The MC simulations are used to correct the distributions at the level of reconstructed particles back to the hadron level on a bin-by-bin basis. The correction factors are determined from the MC simulations as the ratios of the normalised cross sections obtained from particles at hadron level without QED radiation to the normalised cross sections calculated using reconstructed particles and including QED radiation effects. For the forward photon analysis the average of the correction factors determined from LEPTO and CDM is used. For the forward neutron analysis the correction factors are calculated using the combination of RAPGAP-\(\pi \) and CDM simulations, with the weighting factors \(0.6\) and \(1.4\), as described above. The size of the correction factors varies between \(2\) and \(3.5\) for the forward photon and between \(2\) and \(6\) for the forward neutron \(x_F\) distributions, and is about \(3.2\) for the \(W\) distribution in both cases. The correction factors are large mainly due to the non-uniform azimuthal acceptance of the FNC, which is about \(30~\%\) on average. The bin purity, defined as the fraction of events reconstructed in a particular bin that originate from the same bin on hadron level, varies between 75 and \(95~\%\).

2.6 Systematic uncertainties

Several sources of experimental uncertainties are considered and their effect on the measured cross section is quantified. The systematic uncertainties on the cross section measurements are determined using MC simulations, by propagating the corresponding uncertainty through the full analysis chain.

As the cross sections are normalised to the inclusive DIS cross section measured in this analysis, some important systematic uncertainties, such as those involving the trigger efficiency and the integrated luminosity and those related to the reconstruction of the scattered positron and the hadronic final state are largely reduced or cancel. The following sources are considered for both the DIS sample and the FNC samples:

  • The uncertainty on the measurements of the scattered positron energy (\(1~\%\)).

  • The uncertainty on the measurements of the scattered positron angle (\(1\) mrad).

  • The uncertainty on the measurements of the energy of the hadronic final state (\(4~\%\)).

  • The uncertainty on the trigger efficiency (\(1~\%\)).

These uncertainties are strongly correlated between the DIS and the forward photon and neutron samples. The resulting combined uncertainty of the cross section is about \(2~\%\) on average and is considered as uncorrelated between the measurement points.

Several sources of uncertainties related to the reconstruction of the forward photons and neutrons in the FNC are considered:

  • The acceptance of the FNC calorimeter is defined by the interaction point and the geometry of the HERA magnets and is determined using MC simulations. The uncertainty of the impact position of the particle on the FNC, due to beam inclination and the uncertainty on the FNC position, is estimated to be \(5\) mm. This results in uncertainties on the FNC acceptance determination of up to \(15~\%\) for the \(x_F\) distributions.

  • The absolute electromagnetic energy scale of the FNC is known to a precision of \(5~\%\) [8]. This leads to an uncertainty of \(1~\%\) on the forward photon cross section measurement at low energies, increasing to \(11~\%\) for the largest \(x_F\) values.

  • The uncertainty in the neutron detection efficiency and the \(2~\%\) uncertainty on the absolute hadronic energy scale of the FNC [7] lead to systematic errors on the cross section of \(2~\%\) and up to \(10~\%\), respectively.

These systematic uncertainties related to the FNC are strongly correlated between measurement bins and mainly contribute to the overall normalisation uncertainty. For the normalised cross sections studied as a function of \(W\) all above-mentioned FNC related systematic uncertainties contribute to a normalisation uncertainty of approximately \(7~\%\).

In the procedure of correcting the measurements to the hadron level, using MC simulations, the following sources of systematic uncertainties are considered:

  • The systematic uncertainty arising from the radiative corrections and the model dependence of the data correction for the forward neutron cross section is estimated by varying the DJANGOH and RAPGAP-\(\pi \) scaling factors within values permitted by the data and by switching between the CDM and LEPTO models within DJANGOH. The resulting uncertainty on the cross section is typically \(4~\%\), increasing to \(5~\%\) at lowest and highest \(x_F\) values.

  • For the forward photon cross section the systematic uncertainty, taken as the difference of the acceptance corrections calculated using the LEPTO and CDM models, increases from \(1~\%\) at low \(x_F\) to \(7~\%\) at higher \(x_F\).

  • The use of different parton distribution functions in the MC simulation results in a negligible change of the correction factors.

  • A \(2~\%\) uncertainty is attributed to the bin-by-bin subtraction of background arising from the wrong identification of photon and neutron candidates in the FNC and from photoproduction processes.

These uncertainties are assumed to be equally shared between correlated and uncorrelated parts.

The systematic uncertainties shown in the figures and tables are calculated using the quadratic sum of all contributions, which may vary from point to point. They are significantly larger than the statistical uncertainties in all measurement bins. The total systematic uncertainty for the normalised cross section measurements ranges between 8 and \(22~\%\) for the \(x_F\) dependent cross sections and is about \(8~\%\) for the \(W\) dependent cross sections.

3 Results

3.1 Normalised cross sections as a function of \(W\)

The ratios of the forward photon and forward neutron production cross sections to the inclusive DIS cross section, \(\sigma ^{\gamma ,n}_{DIS}/\sigma _{DIS}\), in the kinematic range \(6<Q^2<100~\mathrm GeV^2\) and \(0.05<y<0.6\) and as a function of \(W\) are presented in Tables 2 and 3 and are shown in Fig. 2. Within uncertainties the measured ratios are consistent with constant values of about \(0.027\) (forwards photons) and \(0.083\) (forward neutrons). In other words, within uncertainties the \(W\) dependence of the cross section is independent of the presence of a forward neutron or a forward photon, as predicted by the limiting fragmentation hypothesis [9, 10].

Table 2 The fraction of DIS events with forward photons in the kinematic region given in Table 1. For each measurement, the statistical (\(\delta _{stat.}\)), the total systematic (\(\delta _{total sys.}\)), the uncorrelated (\(\delta _{uncorrel.sys.}\)) systematic uncertainties and the bin-to-bin correlated systematic uncertainties due to the FNC absolute energy scale (\(\delta _{E_{FNC}}\)), the measurement of the particle impact position in the FNC (\(\delta _{XY_{FNC}}\)) and the model dependence of the data correction (\(\delta _{model}\)) are given
Table 3 The fraction of DIS events with forward neutrons in the kinematic region given in Table 1. For each measurement, the statistical (\(\delta _{stat.}\)), the total systematic (\(\delta _{total sys.}\)), the uncorrelated (\(\delta _{uncorrel.sys.}\)) systematic uncertainties and the bin-to-bin correlated systematic uncertainties due to the FNC absolute energy scale (\(\delta _{E_{FNC}}\)), the measurement of the particle impact position in the FNC (\(\delta _{XY_{FNC}}\)) and the model dependence of the data correction (\(\delta _{model}\)) are given
Fig. 2
figure 2

The fraction of DIS events with forward photons (a) and forward neutrons (b) as a function of \(W\) in the kinematic region given in Table 1. Also shown are the predictions of the LEPTO (solid line) and CDM (dashed line) MC models. In the case of forward neutron production, the predictions of RAPGAP-\(\pi \) and the linear combinations of LEPTO and RAPGAP-\(\pi \), as well as CDM and RAPGAP-\(\pi \) are also shown

In Fig. 2 the MC model calculations are compared to the measurements. Both CDM and LEPTO predict a forward photon rate of about \(70~\%\) higher than observed. A similar excess was observed earlier [8]. The photon production rate as a function of \(W\) is rather flat in CDM and shows a slight increase with \(W\) in LEPTO. The shape of the \(W\) distribution is in both models consistent with the data, within errors.

The rate of forward neutron production predicted by LEPTO is consistent with the data, while CDM predicts a much lower rate. However, as was shown in the previous measurement [7], the energy distribution of forward neutrons can be described by MC simulation only if this includes contributions both from standard fragmentation as simulated in DJANGOH, and from a pion exchange mechanism as explicitly simulated in RAPGAP-\(\pi \) but not included in DJANGOH. In Fig. 2b the combinations of the RAPGAP-\(\pi \) and DJANGOH simulations, as described in Sect. 2.5, are compared to the measurement. The weighting factors \(1.4\) for the CDM, \(0.7\) for the LEPTO and \(0.6\) for the RAPGAP-\(\pi \) predictions are determined by fitting the observed neutron energy distributions integrated over the full \(W\) range. The cross sections for inclusive DIS events, used for the normalisation of the forward neutron cross sections, \(\sigma _{DIS}\), are taken from the CDM and LEPTO simulations without additional weights. The model combination describes the observed \(W\) dependence well. It is remarkable that the factors for the CDM and LEPTO contributions differ by a factor two (\(1.4\) and \(0.7\), respectively). It is also notable that the CDM model, which overestimates the rate of forward photons by about \(70~\%\), has to be scaled up in the combination to describe the forward neutron data.

In Fig. 3 predictions of various cosmic ray hadronic interaction models (EPOS LHC, SIBYLL 2.1 and the two versions of QGSJET) are compared to the measured normalised cross sections as a function of \(W\). The CR model predictions show significant differences in absolute values, for both forward photons and forward neutrons. For photons all models predict too high rates by \(30\) to \(40~\%\), and these rates, with the exception of EPOS LHC, show a slight decrease with increasing \(W\), not confirmed by data. For forward neutrons all CR predictions show a \(W\) independent behaviour, in accordance with the measured \(W\) dependence. The QGSJET 01 model predicts a much too high and SIBYLL 2.1 a much too low neutron rate, while the EPOS LHC and QGSJET II-04 models are closer to the measurement.

Fig. 3
figure 3

The fraction of DIS events with forward photons (a) and forward neutrons (b) as a function of \(W\) in the kinematic region given in Table 1. Also shown are the predictions of the cosmic ray hadronic interaction models SIBYLL 2.1 (solid line), QGSJET 01 (dashed line), QGSJET II-04 (dotted line) and EPOS LHC (dash-dotted line)

3.2 Normalised cross sections as a function of \(x_F\) and test of Feynman scaling

The measured normalised differential cross sections, \( {1}/\sigma _{DIS}~ \mathrm{d}\sigma /\mathrm{d}x_F\), of the most energetic photon are presented as a function of \(x_F\) in Table 4 and in Figs. 4 and 5 for the kinematic region defined in Table 1. In order to study the energy dependence of the \(x_F\) distributions, these cross sections are measured in three \(W\) intervals.

Table 4 Normalised cross sections of forward photon production in DIS as a function of \(x_F\). The kinematic phase space of the measurements is given in Table 1. For each measurement, the statistical (\(\delta _{stat.}\)), the total systematic (\(\delta _{total sys.}\)), the uncorrelated (\(\delta _{uncorrel.sys.}\)) systematic uncertainties and the bin-to-bin correlated systematic uncertainties due to the FNC absolute energy scale (\(\delta _{E_{FNC}}\)), the measurement of the particle impact position in the FNC (\(\delta _{XY_{FNC}}\)) and the model dependence of the data correction (\(\delta _{model}\)) are given
Fig. 4
figure 4

Normalised cross sections of forward photon production in DIS as a function of \(x_F\) in three \(W\) intervals in the kinematic region given in Table 1. The error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of the LEPTO (solid line) and CDM (dashed line) MC models

Fig. 5
figure 5

Normalised cross sections of forward photon production in DIS as a function of \(x_F\) in three \(W\) intervals in the kinematic region given in Table 1. The inner error bars show the statistical uncertainty, while the outer error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of the cosmic ray hadronic interaction models SIBYLL 2.1 (solid line), QGSJET 01 (dashed line), QGSJET 01 (no mi) (dash-double dotted line), QGSJET II-04 (dotted line) and EPOS LHC (dash-dotted line). In the right column the ratios of the CR model predictions to the data are shown

The normalised differential cross sections as a function of \(x_F\) are similar for the three \(W\) ranges. As shown in Fig. 4 and already seen in the comparison of the \(W\) dependence, the LEPTO and CDM models predict a rate of forward photons about \(70~\%\) higher than measured. The shapes of the measured distributions are well described by LEPTO, while the CDM description is very poor by showing a significantly harder spectrum than observed in data. In Fig. 5 the predictions of the CR hadronic interaction models are compared to the same measurements. Large differences between the CR models are observed, both in shape and in normalisation. All models tested here overestimate the forward photon rate by 30 to \(40~\%\) at low \(x_F\). The EPOS LHC model describes the shapes of the photon \(x_F\) distributions well. The SIBYLL 2.1 model predicts a harder \(x_F\) dependence, while the spectra obtained from the different variants of QGSJET are softer than observed in the data. Forward photon and neutral pion measurements at the LHC also revealed differences with respect to a similar selection of CR models [6870].

The normalised differential cross sections for forward neutrons are presented in Table 5 and in Figs. 6 and 7 for the kinematic region defined in Table 1. The \(x_F\) distributions are well reproduced by a combination of CDM and RAPGAP-\(\pi \), using the weighting factors and normalisation as described in Sect. 3.1. The individual contributions of the two models are shown in Fig. 6 as well. Fragmentation, as simulated by CDM, dominates the neutron production at lower \(x_F\), while the contribution from pion exchange becomes significant at \(x_F \gtrsim 0.7\). The combination of LEPTO and RAPGAP-\(\pi \) (not shown) also provides a good description of the measurements for the three \(W\) ranges. In Fig. 7 the predictions of the CR hadronic interaction models are compared to the forward neutron production cross sections. The EPOS LHC model provides a reasonable description of the neutron \(x_F\) distributions, except at the highest \(x_F\) values. The SIBYLL 2.1 model describes the shape of the \(x_F\) spectra but fails in the absolute rate. The QGSJET II-04 model shows a harder \(x_F\) dependence, and QGSJET 01 predicts a much too high neutron rate.

Table 5 Normalised cross sections of forward neutron production in DIS as a function of \(x_F\). The kinematic phase space of the measurements is given in Table 1. For each measurement, the statistical (\(\delta _{stat.}\)), the total systematic (\(\delta _{total sys.}\)), the uncorrelated (\(\delta _{uncorrel.sys.}\)) systematic uncertainties and the bin-to-bin correlated systematic uncertainties due to the FNC absolute energy scale (\(\delta _{E_{FNC}}\)), the measurement of the particle impact position in the FNC (\(\delta _{XY_{FNC}}\)) and the model dependence of the data correction (\(\delta _{model}\)) are given
Fig. 6
figure 6

Normalised cross sections of forward neutron production in DIS as a function of \(x_F\) in three \(W\) intervals in the kinematic region given in Table 1. The inner error bars show the statistical uncertainty, while the outer error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of CDM (dotted line), RAPGAP-\(\pi \) (dashed line) and a linear combination of CDM and RAPGAP-\(\pi \) predictions (solid line)

Fig. 7
figure 7

Normalised cross sections of forward neutron production in DIS as a function of \(x_F\) in three \(W\) intervals in the kinematic region given in Table 1. The inner error bars show the statistical uncertainty, while the outer error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of the cosmic ray hadronic interaction models SIBYLL 2.1 (solid line), QGSJET 01 (dashed line), QGSJET 01 (no mi) (dash-double dotted line), QGSJET II-04 (dotted line) and EPOS LHC (dash-dotted line). In the right column the ratios of the CR model predictions to the data are shown

A modified version of the QGSJET 01 model, denoted ’QGSJET 01 (no mi)’ [67], in which the contribution of multi-parton interactions is excluded (see Sect. 2.5), is also compared to the measurements. When multi-parton interactions are switched off, the predicted \(x_F\) spectra become harder without improving the data description.

The \(W\) dependence of the \(x_F\) distributions allows a test of the Feynman scaling hypothesis for particle production. For this test, the ratios of the normalised cross sections for different CM energy intervals are studied as a function of \(x_F\). Figs. 8 and 9 show the ratios of the second to the first and the third to the first \(W\) range for photons. The predictions from CDM, LEPTO and the CR models are also shown. In Fig. 10 the same ratios are shown for forward neutrons and the CR models are compared to the data. For all data distributions the values of these ratios are consistent with unity and with being constant within uncertainties, suggesting that Feynman scaling in the target fragmentation region holds for photons and neutrons. The LEPTO and CDM MC models, used for the comparison with forward photon data, show a similar behaviour. All CR models indicate deviations from scaling for the forward photons, such that the production rate decreases with increasing \(W\). In particular, this effect is strong for the SIBYLL 2.1 and QGSJET 01 models. For forward neutrons the CR models are consistent with Feynman scaling, with exception of SIBYLL 2.1.

Fig. 8
figure 8

Ratios of normalised cross sections of forward photon production in DIS corresponding to two different \(W\) intervals, shown in Fig. 4, as a function of \(x_F\): (a) ratio of the cross section in the \(130<W<190\) GeV interval to the cross section in the \(70<W<130\) GeV interval; (b) ratio of the cross section in the \(190<W<245\) GeV interval to the cross section in the \(70<W<130\) GeV interval. The kinematic phase space is defined in Table 1. The error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of the LEPTO (solid line) and CDM (dashed line) MC models

Fig. 9
figure 9

Ratios of normalised cross sections of forward photon production in DIS corresponding to two different \(W\) intervals, shown in Fig. 5, as a function of \(x_F\): a ratio of the cross section in the \(130<W<190\) GeV interval to the cross section in the \(70<W<130\) GeV interval; b ratio of the cross section in the \(190<W<245\) GeV interval to the cross section in the \(70<W<130\) GeV interval. The kinematic phase space is defined in Table 1. The error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of the cosmic ray hadronic interaction models SIBYLL 2.1 (solid line), QGSJET 01 (dashed line), QGSJET II-04 (dotted line) and EPOS LHC (dash-dotted line)

Fig. 10
figure 10

Ratios of normalised cross sections of forward neutron production in DIS corresponding to two different \(W\) intervals, shown in Fig. 7, as a function of \(x_F\): a ratio of the cross section in the \(130<W<190\) GeV interval to the cross section in the \(70<W<130\) GeV interval; b ratio of the cross section in the \(190<W<245\) GeV interval to the cross section in the \(70<W<130\) GeV interval. The kinematic phase space is defined in Table 1. The error bars show the total experimental uncertainty, calculated using the quadratic sum of the statistical and systematic uncertainties. Also shown are the predictions of the cosmic ray hadronic interaction models SIBYLL 2.1 (solid line), QGSJET 01 (dashed line), QGSJET II-04 (dotted line) and EPOS LHC (dash-dotted line)

4 Summary

The production of high energy forward neutrons and photons has been studied at HERA in deep-inelastic \(ep\) scattering in the kinematic region \(6< Q^2 < 100~\mathrm {GeV}^2\) and \(0.05<y<0.6\). The normalised DIS cross sections \(1/\sigma _{DIS}~\mathrm{d}\sigma /\mathrm{d}x_F\) for the production of photons and neutrons at pseudorapidities \(\eta >7.9\) and in the range of Feynman-\(x\) of \(0.1<x_F<0.7\) for the photons and \(0.1<x_F<0.94\) for neutrons are presented. The measured cross sections as a function of \(x_F\) at different centre-of-mass energies of the virtual photon-proton system agree within uncertainties, confirming the validity of Feynman scaling in the energy range of the virtual photon-proton system \(70<W<245~\mathrm {GeV}\).

Different Monte Carlo models are compared to the measurements. All these models overestimate the rate of photons by 30–70 %. The shapes of the measured forward photon cross sections are well described by the LEPTO MC simulation, while predictions based on the colour dipole model fail, especially at high \(x_F\). The cross sections for forward neutrons are well described by a linear combination of the standard fragmentation model, as implemented in DJANGOH, and the one-pion-exchange model RAPGAP-\(\pi \). Predictions of models, which are commonly used for the simulation of cosmic ray cascades, are also compared to the forward photon and neutron measurements. None of the models describes the photon and neutron data simultaneously well. The best description of the shapes of the photon and the neutron \(x_F\) distributions is provided by the EPOS LHC model. Within the kinematic range of the measurements, the relative rate of forward photons and neutrons in DIS events is observed to be independent of the energy of the virtual photon-proton CM, and therefore also consistent with the hypothesis of limiting fragmentation. The present measurement provides new information to further improve the understanding of proton fragmentation in collider and cosmic ray experiments.