1 Introduction

Measurements of leading baryon production in high energy particle collisions, i.e. the production of protons and neutrons at very small polar angles with respect to the initial hadron beam direction (forward direction), are important inputs for the theoretical understanding of strong interactions in the soft, non-perturbative regime. In ep collisions at HERA, a hard scale may be present in such reactions if the photon virtuality, \(Q^2\), is large, or if objects with high transverse momenta, \(p_T\), are produced in addition to the leading baryon. In such cases the process usually can be factorised into short-distance and long-distance phenomena and perturbative QCD often is applicable for the description of the hard part of the process.

Previous HERA measurements [17] have demonstrated that in the semi-inclusive reaction \( e+p \rightarrow e+n+X\) the production of neutrons carrying a large fraction of the proton beam energy is dominated by the pion exchange process. In this picture a virtual photon, emitted from the beam electron, interacts with a pion from the proton cloud, thus giving access to the \(\gamma ^*\pi \) cross section and, in the deep-inelastic scattering regime, to the pion structure function.

The aim of the present analysis is to measure exclusive \(\rho ^0\) production on virtual pions in the photoproduction regime at HERA and to extract the quasi-elastic \(\gamma \pi \rightarrow \rho ^0\pi \) cross section for the first time. Since no hard scale is present, a phenomenological approach, such as Regge theory [8], is most appropriate to describe the reaction. In the Regge framework such events are explained by the diagram shown in Fig. 1a which involves an exchange of two Regge trajectories in the process \(2 \rightarrow 3\), known as a Double Peripheral Process (DPP), or Double-Regge-pole exchange reaction [911]. This process can also be seen as a proton dissociating into \((n,\pi ^+)\) system which scatters elastically on the \(\rho ^0\) via the exchange of the Regge trajectory with the vacuum quantum numbers, called the “Pomeron”.

In the past, similar reactions were studied at lower energies in nucleon–nucleon and meson–nucleon collisions [1217]. Most of the experimental properties of these reactions were successfully explained by the generalised Drell–Hiida–Deck model (DHD) [1821], in which in addition to the pion exchange (Fig. 1a) two further contributions (Fig. 1b, c) are included. The graphs depicted in Fig. 1b, c give contributions to the total scattering amplitude with similar magnitude but opposite sign [22, 23]. Therefore they largely cancel in most of the phase space, in particular at small momentum transfer squared at the proton vertex, \(t \rightarrow 0\), such that the pion exchange diagram dominates the cross section [21]. One of the specific features observed in these experiments is a characteristic \(t^{\prime }\) dependence at the ‘elastic’ vertex,Footnote 1 with the slope dependent on the mass of the \((n\pi )\) system produced at the other, \(pn\pi ^+\), vertex, and changing in a wide range of approximately \(4<b(m)<22\) GeV\(^{-2}\). The Deck model in its original formulation cannot fully describe such a strong mass-slope correlation and interference between the amplitudes corresponding to the first three graphs in Fig. 1 has to be taken into account to explain the experimental data [24, 25].

In the analysis presented here only the two charged pions from the \(\rho ^0\) decay and the leading neutron are observed directly. The pion from the proton vertex is emitted under very small angles with respect to the proton beam and escapes detection. This leads to a background contamination from events with a different final state, which originate from diffractive dissociation of the proton into a system Y containing a neutron (Fig. 1d). Using the H1 detector capabilities in the forward region such processes can be suppressed to a certain extent. The residual background contribution is estimated from a Monte Carlo model tuned to describe vector meson production in diffractive dissociation at HERA.

Fig. 1
figure 1

Generic diagrams for processes contributing to exclusive photoproduction of \(\rho ^0\) mesons associated with leading neutrons at HERA. The signal corresponds to the Drell–Hiida–Deck model graphs for the pion exchange (a), neutron exchange (b) and direct pole (c). Diffractive scattering in which a neutron may be produced as a part of the proton dissociation system, \(M_Y\), contributes as background (d). The \(N^*\) in c denotes both resonant (via \(N^+\)) and possible non-resonant \(n+\pi ^+\) production

The analysis is based on a data sample corresponding to an integrated luminosity of 1.16 pb\(^{-1}\) collected with the H1 detector in the years 2006 and 2007. 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}\). The photon virtuality is limited to \(Q^2<2~\mathrm {GeV}^2\) with an average value of \(0.04~\mathrm {GeV}^2\).

2 Cross sections definitions

The kinematics of the process

$$\begin{aligned} e(k) + p(P) \rightarrow e(k') + \rho ^0(V) + n(N) + \pi ^+, \end{aligned}$$
(1)

where the symbols in parentheses denote the four-momenta of the corresponding particles, is described by the following invariants:

  • the square of the ep centre-of-mass energy \(s = (P+k)^2\),

  • the modulus of the four-momentum transfer squared at the lepton vertex \(Q^2 = -q^2 = -(k-k')^2\),

  • the inelasticity \(y=(q \cdot P)/(k \cdot P)\),

  • the square of the \(\gamma p\) centre-of-mass energy \(W_{\gamma p}^2 = (q+P)^2 \simeq ys-Q^2\),

  • the fraction of the incoming proton beam energy carried by the leading neutron \(x_L = (q \cdot N)/(q \cdot P) \simeq E_n/E_p\),

  • the four-momentum transfer squared at the proton vertex \(t = (P-N)^2 \simeq -\frac{p_{T,n}^2}{x_L} - \frac{(1-x_L)(m_n^2-m_p^2x_L)}{x_L}\), and

  • the four-momentum transfer squared at the photon vertex \(t' = (q-V)^2\).

Here \(E_p, m_p, E_n, m_n\) represent the energy and the mass of the incoming proton and the outgoing leading neutron, respectively. The square of the \(\gamma \pi \) centre-of-mass energy is then given by \(W_{\gamma \pi }^2 \simeq W_{\gamma p}^2(1-x_L)\).

Experimentally, the kinematic variables at the photon vertex (the mass \(M_{\rho }\), the pseudorapidity \(\eta _{\rho }\) and the transverse momentum squared \(p_{T,\rho }^2\) of the \(\rho ^0\) meson) are determined from the \(\rho ^0\) decay pions, while those at the proton vertex (\(x_L\) and \(p_{T,n}^2\)) are deduced from the measured energy and scattering angle of the leading neutron.

In the limit of photoproduction, i.e. \(Q^2 \rightarrow 0\), the beam positron is scattered at small angles and escapes detection. In this regime the square of the \(\gamma p\) centre-of-mass energy can be reconstructed via the variable \(W_{\gamma p,rec}^2 = s \, y^{rec}\), where \(y^{rec}\) is the reconstructed inelasticity, measured as \(y^{rec} = (E_{\rho }-p_{z,\rho })/(2E_e)\). Here, \(E_{\rho }\) and \(p_{z,\rho }\) denote the reconstructed energy and the momentum along the proton beam direction (z-axis) of the \(\rho ^0\) meson and \(E_e\) is the positron beam energy. The variable \(t'\) can be estimated from the transverse momentum of the \(\rho ^0\) meson in the laboratory frame via the observable \(t'_{rec} = -p_{T,\rho }^2\) to a very good approximation.Footnote 2

The cross section of the exclusive reaction (1) can be expressed as a product of a virtual photon flux \(f_{\gamma /e}\) and a photon–proton cross section \(\sigma _{\gamma p}\):

$$\begin{aligned} \frac{\mathrm{d}^2\sigma _{ep}}{\mathrm{d}y \mathrm{d}Q^2} = f_{\gamma /e}(y,Q^2) \sigma _{\gamma p}(W_{\gamma p}(y)). \end{aligned}$$
(2)

In the Vector Dominance model (VDM) [2628] taking into account both transversely and longitudinally polarised virtual photons the effective photon flux is given by

$$\begin{aligned}&f_{\gamma /e}(y,Q^2) = \frac{\alpha }{2\pi Q^2 y}\nonumber \\&\quad \times \left\{ \left[ 1+(1-y)^2 - 2(1-y) \left( \frac{Q^2_\mathrm{min}}{Q^2}-\frac{Q^2}{M^2_{\rho }} \right) \right] \frac{1}{\left( 1+\frac{Q^2}{M^2_{\rho }}\right) ^2} \right\} , \end{aligned}$$
(3)

where \(\alpha \) is the fine structure constant and \(Q_{\min }^2=m_e^2 y^2/(1-y)\), with \(m_e\) being the mass of the electron and \(M_{\rho }\) is the \(\rho ^0\) meson mass.

In the one-pion-exchange (OPE) approximation [2931], which is valid for very small \(p_{T,n}^2 \sim m_{\pi }^2\), the photon–proton cross section can be further decomposed into a pion flux, describing \(p \rightarrow n\pi ^+\) splitting, convoluted with a photon-pion cross section:

$$\begin{aligned} \frac{\mathrm{d}^2\sigma _{\gamma p}(W_{\gamma p},x_L,t)}{\mathrm{d}x_L\mathrm{d}t} = f_{\pi /p}(x_L,t)\, \sigma _{\gamma \pi }(W_{\gamma \pi }). \end{aligned}$$
(4)

A generic expression for the pion flux factor can be written as follows:

$$\begin{aligned} f_{\pi /p}(x_L,t)= & {} \frac{1}{2\pi } \frac{g_{p\pi n}^2}{4\pi } (1-x_L)^{\alpha _{I\!\!P}(0)-2\alpha _{\pi }(t)}\nonumber \\&\times \frac{-t}{(m_{\pi }^2-t)^2} F^2(t,x_L), \end{aligned}$$
(5)

where \(\alpha _{I\!\!P}(0)\) is the Pomeron intercept, \(\alpha _{\pi }(t)=\alpha ^{\prime }_{\pi }(t-m_{\pi }^2)\) is the pion trajectory, \(g_{p\pi n}^2/4\pi \) is the \(p\pi n\) coupling constant known from phenomenological analysis [32] of low energy data, and \(F(t,x_L)\) is a form factor accounting for off mass-shell corrections and normalised to unity at the pion pole, \(F(m_{\pi }^2,x_L)=1\). There exists a variety of models for the exact form of the pion flux [3340] which typically leads to a \({\sim }30~\%\) spread in the predicted cross section according to Eq. (4). Most of models use a non-Reggeized version of Eq. (5), i.e. \(\alpha _{I\!\!P}(0)=1\) and \(\alpha _{\pi }(t)=0\).

3 Experimental procedure and data analysis

3.1 H1 detector

A detailed description of the H1 detector can be found elsewhere [41,42]. Only those components relevant for the present analysis are described here. The origin of the right-handed H1 coordinate system is the nominal ep 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 xy plane. The pseudorapidity is defined by \(\eta = -\ln {[\tan (\theta /2)]}\) and is measured in the laboratory frame.

The central region of the detector is equipped with a tracking system. It included a set of two large coaxial cylindrical drift chambers (CJC), interleaved by a z chamber, and the central silicon tracker (CST) [42] operated in a solenoidal magnetic field of \(1.16~\mathrm{T}\). This provides a measurement of the transverse momentum of charged particles with resolution \(\sigma (p_T) / p_T \simeq 0.002 \ p_T \oplus 0.015\) (\(p_T\) measured in GeV), for particles emitted from the nominal interaction point with polar angle \( 20^{\circ } \le \theta \le 160^{\circ }\). The interaction vertex is reconstructed from the tracks. The five central inner proportional chambers (CIP) [43] are located between the inner CJC and the CST. The CIP has an angular acceptance in the range \(10^{\circ }<\theta <170^{\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 is used to measure the scattered electron and to reconstruct the energy of the hadronic final state. The backward region (\(153^\circ <\theta <177.8^\circ \)) is covered by a lead/scintillating-fibre calorimeter (SpaCal) [44]; its main purpose is the detection of scattered positrons.

A set of “forward detectors” is sensitive to the energy flow close to the outgoing proton beam direction. It consists of the forward muon detector (FMD), the Plug calorimeter and the forward tagging system (FTS). The lead–scintillator Plug calorimeter enables energy measurements to be made in the pseudorapidity range \(3.5 < \eta < 5.5\). It is positioned around the beam-pipe at \(z = 4.9\) m. The FMD is a system of six drift chambers which are grouped into two three-layer sections separated by a toroidal magnet. Although the nominal coverage of the FMD is \(1.9 < \eta < 3.7\), particles with pseudorapidity up to \(\eta \simeq 6.5\) can be detected indirectly through their interactions with the beam transport system and detector support structures. The very forward region, \(6.0 < \eta < 7.5\), is covered by an FTS station which is used in this analysis. It consists of scintillator detectors surrounding the beam pipe at \(z=28\) m. The forward detectors together with the LAr calorimeter are used here to suppress inelastic and proton dissociative background by requiring a large rapidity gap (LRG) void of activity between the leading neutron and the pions from the \(\rho ^0\) decay.

Neutral particles produced at very small polar angles can be detected in the forward neutron calorimeter (FNC) [7, 45], which is situated at 106 m from the interaction point. It covers the pseudorapidity range \(\eta > 7.9\). 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 hadronic interaction lengths \(\lambda \), and the Main Calorimeter with a total length of \(8.9 \lambda \) (see Fig. 2a). 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, as illustrated in Fig. 2b.

Fig. 2
figure 2

Sketch of the Forward Neutron Calorimeter (a), acceptance in the azimuthal plane (b) and the \((x_L,p_{T,n}^2)\) plane (c). The shaded area in b is the projected aperture limited by the proton beamline elements. The inset of c shows the acceptance in terms of the \(x_L\) and t variables. The dark green area indicates the OPE2 region, \(p_{T,n}<0.2\) GeV and \(x_L>0.65\). The curve in the main figure corresponds to the angular cut \(\theta _n<0.75\) mrad, and the grid shows the binning scheme used for the double differential cross section measurement, \(\mathrm{d^2}\sigma _{\gamma p}/\mathrm{d}x_L\mathrm{d}p_{T,n}^2\). The dots are events from the preselection sample described in the leftmost column of Table 1

The absolute electromagnetic and hadronic energy scales of the FNC are known to 5 and \(2~\%\) precision, respectively [7]. 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. The spatial resolution is \(\sigma (x,y)\approx 10~\hbox {cm}/\sqrt{E~[\mathrm GeV]} \oplus 0.6~\hbox {cm}\) for hadronic showers starting in the Main Calorimeter. A better spatial resolution of about \(2~\hbox {mm}\) is achieved for electromagnetic showers and for those hadronic showers which start in the Preshower Calorimeter.

The instantaneous luminosity is monitored based on the rate of the Bethe–Heitler process \(ep \rightarrow ep\gamma \). The final state photon is detected in the photon detector located close to the beampipe at \(z=-103\) m. The precision of the integrated luminosity measurement is improved in a dedicated analysis of the elastic QED Compton process [46] in which both the scattered electron and the photon are detected in the SpaCal.

3.2 Event selection

The data sample of this analysis has been collected using a special low multiplicity trigger requiring two tracks with \(p_T>160\) MeV and originating from the nominal event vertex, and at most one extra track with \(p_T>100\) MeV. The tracks are found by the Fast Track Trigger (FTT) [47, 48], based on hit information provided by the CJCs. The trigger also contains a veto condition against non-ep background provided by the CIP. The average trigger efficiency is about \(75 ~\%\) for the analysis phase space. The trigger simulation has been verified and tuned to the data using an independently triggered data sample.

For the analysis, exclusive events are selected, containing two oppositely charged pion candidates in the central tracker, a leading neutron in the FNC and nothing else above noise level in the detector.Footnote 3 The photoproduction regime is ensured by the absence of a high energy electromagnetic cluster consistent with a signal from a scattered beam positron in the calorimeters. This limits the photon virtuality to \(Q^2 \lesssim 2 \, \mathrm {GeV}^2 \), resulting in a mean value of \(\langle Q^2 \rangle = 0.04 \, \mathrm {GeV}^2 \).

The \(\rho ^0\) candidate selection requires the reconstruction of the trajectories of two, and only two, oppositely charged particles in the central tracking detector. They must originate from a common vertex lying within \(\pm 30\) cm in z of the nominal ep interaction point, and must have transverse momenta above 0.2 GeV and polar angles within the interval \( 20^{\circ } \le \theta \le 160^{\circ }\). The momentum of the \(\rho ^0\) meson is calculated as the vector sum of the two charged particle momenta. The two-pion invariant mass is required to be within the interval \(0.3 < M_{\pi \pi }<1.5\) GeV. Since no explicit hadron identification is used, events are discarded with \(M_{KK}<1.04\) GeV where \(M_{KK}\) is the invariant mass of two particles under the kaon mass hypothesis. This cut suppresses a possible background from exclusive production of \(\phi \) mesons.

Events containing a leading neutron are selected by requiring a hadronic cluster in the FNC with an energy above \(120~\mathrm {GeV}\) and a polar angle below \(0.75~\mathrm {mrad}\). The cut on polar angle, defined by the geometrical acceptance of the FNC, restricts the neutron transverse momenta to the range \(p_{T,n}< x_L \cdot 0.69~\mathrm {GeV}\).

To ensure exclusivity, additional cuts are applied on the calorimetric energy and on the response of the forward detectors. There should be no cluster with energy above 400 MeV, unless associated with \(\rho ^0\) decay products, in the SpaCal and LAr calorimeters. A Large Rapidity Gap signature is required, by selecting events with no activity above noise levels in the forward detectors. This suppresses non-diffractive interactions to a negligible level and also significantly reduces diffractive background.

Table 1 Event selection criteria and the definition of the kinematic phase space (PS) of the measurements. The measured cross sections are determined at \(Q^2=0\) using the effective flux (3), based on the VDM

After these cuts the data sample contains about 7000 events. The event selection criteria together with the analysis and the measurement phase space definitions are summarised in Table 1. In order to better control migration effects and backgrounds most of the selection cuts are kept softer than the final measurement phase space limits. In the end, the \(\gamma p\) cross sections measured in the \(\theta _n<0.75\) mrad range are based on \({\sim }5770\) events. For the \(\gamma \pi \) cross section extraction additional cuts are applied in order to stay within a range where the validity of OPE can be safely expected. Two sub-samples are defined: OPE1 with \(p_{T,n}<200\) MeV, containing \({\sim }3600\) events and OPE2 with \(p_{T,n}<200\) MeV and \(x_L>0.65\), containing \({\sim }2200\) events.Footnote 4

3.3 Monte Carlo simulations and corrections to the data

Monte Carlo (MC) simulations are used to calculate acceptances and efficiencies for triggering, track reconstruction, event selection and background contributions and to account for migrations between measurement bins due to the finite detector resolution.

Signal events from the DPP reaction (Fig. 1a) are modelled by the two-step MC generator POMPYT [49], in which the virtual pion is produced at the proton vertex according to one of the available pion flux parametrisations. This pion then scatters elastically on the photon from the electron beam, thus producing a vector meson (\(\rho ^0\) in our case). In this analysis the non-Reggeized pion flux factor is taken from the light-cone representation [50] with the form factor in Eq. (5)

$$\begin{aligned} F^2(t,x_L) = \exp \left( - R^2_{\pi n} \frac{m^2_\pi -t}{1-x_L}\right) , \end{aligned}$$
(6)

where \(R_{\pi n}=0.93~\mathrm {GeV}^{-1}\) is the radius of the pion-proton Fock state [34, 35]. The same version of the pion flux factor has been used in previous H1 publications on leading neutron measurements [4, 7] providing a good description of inclusive neutron spectra. For the numerical value of the \(p\pi n\) coupling constant, the most recent estimate [32] \(g_{p\pi n}^2/4\pi = 14.11 \pm 0.20\) is used.

Since the exact shape of the \(p_{T,\rho }^2\) dependence is not a priori known, two extreme versions are generated. In the first version a simple exponential shape is assumed, as expected for elastic \(\rho ^0\) photoproduction on the pion, with the slope \(b=5\) GeV\(^{-2}\). For the second version a mass-dependent slope is taken, \(4\le b(M_{n\pi })\le 22\) GeV\(^{-2}\), typical for DPP processes as observed at lower energies [1214, 21]. The difference in the correction factors obtained using these two versions of MC simulations is part of the model dependent systematic uncertainty.

The background events originating from diffractive \(\rho ^0\) production (Fig. 1d) are generated using the program DIFFVM [51], which is based on Regge theory and the Vector Dominance Model. All channels (elastic, single- and double-dissociation processes) are included, with the relative composition as measured in [52]. For the proton dissociative case the \(M_Y\) mass spectrum is parametrised as \(\mathrm {d} \sigma / \mathrm {d}M_Y^2 \propto 1/M_Y^{2.16}\), for \(M_Y^2 > 3.6~\mathrm {GeV}^2 \) with quark and diquark fragmentation using the JETSET program [53]. For the low mass dissociation the production of excited nucleon states at the proton vertex is taken into account explicitly. Signal events, corresponding to the diagram shown in Fig. 1c, are excluded from the generated background sample.

The DIFFVM program is also used to estimate possible contaminations from diffractive \(\omega (782),~ \phi (1020)\) and \(\rho ^{\prime }(1450{-}1700)\) production.

As discussed in Sect. 1, the pion exchange diagram dominates the cross section in the low t region where the contributions from the diagrams in Fig. 1b, c almost cancel. To check a possible influence of these terms on the MC correction factors, neutron exchange events (b) were generated using POMPYT and events of class (c) using DIFFVM. As expected, these events have kinematic distributions and selection efficiencies similar to those from the pion exchange process and do not alter the MC correction factors beyond the quoted systematic uncertainties.

In both the POMPYT and the DIFFVM generators a simple non-relativistic Breit–Wigner shape is used for the \(\rho \) meson mass. Therefore all MC events are reweighted to the relativistic Breit–Wigner shape with additional \(p_T\)-dependent distortion as observed in \(\rho ^0\) photoproduction experiments. The distortion is caused by the interference between the resonant and non-resonant \(\pi ^+\pi ^-\) production and is characterised by the phenomenological skewing parameter, \(n_{RS}\), as suggested by Ross and Stodolsky [54]:

$$\begin{aligned} \frac{dN(M_{\pi \pi })}{dM_{\pi \pi }} \propto BW_{\rho }(M_{\pi \pi }) \left( \frac{M_{\rho }}{M_{\pi \pi }}\right) ^{n_{RS}(p_{T,\rho })} \end{aligned}$$
(7)

with \(M_{\rho }\) being the nominal resonance mass [55] and \(n_{RS}(p_{T,\rho })\) taken from published ZEUS data on elastic photoproduction of \(\rho ^0\) mesons [56]. Additionally, the signal MC events (POMPYT) are reweighted in \(W_{\gamma p}\) and in \(p_{T,\rho }^2\) to the observed shapes of the corresponding distributions. This reweighting is performed iteratively and has converged after two iterations. The uncertainty in the reweighting procedure is then taken into account in the systematic error analysis.

Small, but non-zero values of \(Q^2\) cause \(|t'|\) to differ from \(p_{T,\rho }^2\) by less than \(Q^2\). To account for this effect a multiplicative correction factor determined with the Monte Carlo generators is applied to the bins of the \(p_{T,\rho }^2\) distribution; the correction is obtained by taking the ratio between the \(|t'|\) and \(p_{T,\rho }^2\) distributions at the generator level. This correction varies from 1.1 at \(p_{T,\rho }^2 = 0\) to 0.77 at \(p_{T,\rho }^2 = 1\) GeV\(^2\).

For all MC samples detector effects are simulated in detail with the GEANT program [57]. The MC description of the detector response, including trigger efficiencies, is adjusted using comparisons with independent data. Beam-induced backgrounds are taken into account by overlaying the simulated events with randomly triggered real events. The simulated MC events are passed through the same reconstruction and analysis chain as is used for the data.

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 size of the correction factors is 12 in average, corresponding to an efficiency of \({\sim }8~\%\), and varies between \({\sim }10\) and \({\sim }24\) for different parts of the covered phase space. The main contributions to the inefficiency are: the azimuthal acceptance of the FNC (\({\sim }30~\%\) on average), the \(\rho \) meson reconstruction efficiency which is zero if one of the tracks has low transverse momentum (\({\sim }60~\%\)), the LRG selection efficiency (\({\sim }60~\%\)) and the trigger efficiency (\({\sim }75~\%\)). 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 70 and \(95~\%\) for one-dimensional distributions and between 45 and \(65~\%\) for two-dimensional ones. As an example, Fig. 2c illustrates the binning scheme used in the two-dimensional \((x_L,p_{T,n})\) distribution.

3.4 Extraction of the \(\rho ^0\) signal

The invariant mass distribution of the two tracks under the charged pion mass hypothesis is shown in Fig. 3a. The distribution is corrected for the mass dependent detector efficiency.

Fig. 3
figure 3

The \(\rho ^0\) meson properties: a mass distribution of the \(\pi ^+\pi ^-\) system for exclusive \(\rho ^0\) production with \(p_T^2<1.0\) GeV\(^2\) associated with a leading neutron. The data points are corrected for the detector efficiency. The curves represent different components contributing to the measured distribution and the Breit–Wigner resonant part extracted from the fit to the data. The analysis region \(0.6<M_{\pi ^+\pi ^-}<1.1\) GeV is indicated by vertical arrows. b Ross–Stodolsky skewing parameter, \(n_{RS}\), as a function of \(p_T^2\) of the \(\pi ^+\pi ^-\) system. The values measured in this analysis are compared to previously obtained results for elastic photoproduction of \(\rho ^0\) mesons, \(\gamma p \rightarrow \rho ^0 p\), by the ZEUS Collaboration. c Decay angular distribution of the \(\pi ^+\) in the helicity frame. d Spin-density matrix element, \(r_{00}^{04}\), as a function of \(Q^2\) for diffractive \(\rho ^0\) photo- and electro-production. The curves on bd represent the results of the fits discussed in the text

A fit is performed in the range \(M_{\pi \pi }>0.4\) GeV using the Ross–Stodolsky parametrisation (7) for the \(\rho ^0\) meson mass shape and adding the contributions for the reflection from \(\omega \rightarrow \pi ^+\pi ^-\pi ^0\) and for the non-resonant background. Other sources of non-\(\rho ^0\) background, such as \(\omega (782) \rightarrow \pi ^+\pi ^-, ~~ \phi (1020)\rightarrow K_L^0 K_S^0, \, \pi ^+\pi ^-\pi ^0, ~~ \rho ^{\prime } \rightarrow \rho \pi \pi , 4\pi , \pi \pi \), which may be misidentified as \(\rho ^0\) candidates, are estimated using MC simulations with the relative yield normalisation fixed to previously measured and published values: \(\sigma _{\gamma p}(\omega )/\sigma _{\gamma p}(\rho ^0)=0.10(\pm 20~\%)\) [58], \(\sigma _{\gamma p}(\phi )/\sigma _{\gamma p}(\rho ^0)=0.07(\pm 20~\%)\) [59] and \(\sigma _{\gamma p}(\rho ^{\prime })/\sigma _{\gamma p}(\rho ^0)=0.20(\pm 50~\%)\) [6063]. The resulting overall background contamination in the analysis region \(0.6<M_{\pi \pi }<1.1\) GeV is found to be \((1.5 \pm 0.7)~\%\).

The fitted values of the resonance mass and width are \(764 \pm 3 (\mathrm stat.)\) MeV and \(155 \pm 5 (\mathrm stat.)\) MeV, respectively, in agreement with the nominal PDG values of \(M_{\rho }\) and \(\Gamma _{\rho }\) [55]. The cross section is then calculated for the full mass range \(2m_{\pi }<M_{\pi \pi }<M_{\rho }+5\Gamma _{\rho }\) using the resonant part only, represented by the relativistic Breit–Wigner function \(BW_{\rho }(M_{\pi \pi })\) with momentum dependent width \(\Gamma (M_{\pi \pi })\) [64]:

$$\begin{aligned}&BW_{\rho }(M_{\pi \pi }) = \frac{M_{\pi \pi } \, M_{\rho } \, \Gamma (M_{\pi \pi })}{(M_{\rho }^2-M_{\pi \pi }^2)^2 + M_{\rho }^2 \, \Gamma (M_{\pi \pi })^2},\nonumber \\&\quad \Gamma (M_{\pi \pi }) = \Gamma _{\rho } \, \left( \frac{q^*}{q_0^*} \right) ^3 \, \frac{M_{\rho }}{M_{\pi \pi }} \end{aligned}$$
(8)

where \(q^*\) is the momentum of the decay pions in the rest frame of a pair of pions with mass \(M_{\pi \pi }\), and \(q^*_0\) is the value of \(q^*\) for \(M_{\pi \pi }=M_{\rho }\).

The Breit–Wigner shape is strongly distorted due to interference with the non-resonant \(\pi \pi \) production amplitude (dashed curve in Fig. 3a). The strength of the distortion is \(p_T\)-dependent and within the ansatz (7) is characterised by the phenomenological skewing parameter, \(n_{RS}\). For the full \(p_T\) range of the present analysis, \(p_{T,\rho }^2<1\) GeV\(^2\), a fit results in the value \(n_{RS} = 4.22 \pm 0.28\). To study its \(p_T\) dependence the fit is repeated in four \(p_{T,\rho }^2\) bins. The values obtained are shown in Fig. 3b in comparison with previously published ZEUS results [56] from elastic \(\rho ^0\) photoproduction, \(\gamma p \rightarrow \rho ^0p\). The dashed curve represents a fit to all these data by the empirical formula

$$\begin{aligned} n_{RS} = n_0 \, (p_T^2 + M^2)^{-\beta } \end{aligned}$$
(9)

with \(n_0, M\) and \(\beta \) as free parameters. The fitted value of \(M^2 \simeq 0.6\) GeV\(^2\) suggests that the relevant scale for photoproduction of vector mesons is indeed \((p_T^2 + M_V^2)\).

An important set of observables which characterise the helicity structure of the vector meson production are the angular distributions of the decay pions. Here we study the distribution of \(\theta _h\) which gives access to the \(\rho ^0\) spin-density matrix element \(r_{00}^{04}\). The angle \(\theta _h\) is defined as the polar angle of the positively charged decay pion in the \(\rho ^0\) rest frame with respect to the meson direction in the \(\gamma ^*p\) centre-of-mass frame. According to the formalism presented in [65] the distribution \(\theta _h\) is given by:

$$\begin{aligned} \frac{1}{\sigma }\frac{\mathrm{d}\sigma }{\mathrm{d}\cos \theta _h} \propto 1-r_{00}^{04} +(3 r_{00}^{04}-1) \cos ^2 \theta _h. \end{aligned}$$
(10)

Figure 3c shows the acceptance corrected \(\cos \theta _h\) distribution together with the fit by Eq. (10) yielding the value of \(r_{00}^{04} = 0.108 \pm 0.017\). In Fig. 3d this result is compared to the values obtained in diffractive \(\rho ^0\) photo- and electro-production at HERA [56, 6669]. The steep \(Q^2\) dependence is driven mainly by the QED gauge invariance motivated factor \(Q^2/M_{V}^2\), and can be fitted by a simple expression [70]

$$\begin{aligned} r_{00}^{04} = \frac{1}{1 + \xi (M_{\rho }^2/Q^2)^{\kappa }} \end{aligned}$$

with the parameters \(\xi = 1.85 \pm 0.10\) and \(\kappa = 0.67 \pm 0.03\), as illustrated by the dashed curve.

In summary, all properties of the selected \(\pi ^+\pi ^-\) sample investigated here are consistent with \(\rho ^0\) photoproduction.

3.5 Signal and background decomposition

The event selection described in Sect. 3.2 does not completely suppress non-DPP background. According to the MC simulations, the remaining part is mostly due to proton dissociation with some admixture of double dissociative events.

As in the case of inclusive leading neutron production [7] signal and background events have different shapes of the leading neutron energy distribution, although in the present analysis the difference is less pronounced. The shape differences in the neutron energy spectrum predicted by MC for the DPP events (POMPYT) and for the proton dissociative background (DIFFVM) are still sufficient to disentangle these two contributions on a statistical basis. For this purpose a combination of the spectra obtained for reconstructed events of these two MC models fulfilling all selection criteria with free normalisation is fitted to the data. From this fit the background fraction is determined to be \(F_{bg}=0.34\pm 0.05\). The uncertainty includes both the fit error and systematic uncertainties related to the background shape variation in terms of \(M_Y\) and t dependencies and proton dissociation fraction in the overall diffractive cross section. Figure 4 illustrates this decomposition using the nominal DIFFVM parameters.

Fig. 4
figure 4

Signal to background decomposition in the selected data sample. On the left panel the distribution of the measured neutron energy, \(E_n\), is shown together with the contributions from signal and background. On the right panel the \(\chi ^2\) dependence on the background fraction, \(F_{bg}\), is shown. The shaded band represents the \(1 \sigma \) uncertainty around the optimal fit value of the \(F_{bg}\), taking into account statistical errors, FNC calibration systematics and the uncertainty in proton dissociation background shape

Fig. 5
figure 5

Distributions of the reconstructed quantities \(p_T^2\) and \(\eta \) of the \(\rho ^0\) meson, \(x_L\), \(p_T\) and \(\varphi \) of the neutron and \(W_{\gamma p}\) for data and Monte Carlo simulations normalised to the data. Data points are shown with statistical errors only. The shaded band indicates the uncertainty in the estimated background fraction

Control plots for the data description by the Monte Carlo models using this signal to background ratio are shown in Fig. 5. Since neither POMPYT nor DIFFVM are able to provide reliable absolute cross section predictions for such a final state, only a shape comparison is possible. The irregular shape of the azimuthal angle distribution, \(\varphi _n\), is due to the FNC aperture limitations, as shown in Fig. 2b.

In the fit described above the absolute normalisation for the DIFFVM prediction is left free. As a cross check, this normalisation has been fixed using an orthogonal, background dominated sample, obtained by requiring an ‘anti-LRG’ selection, i.e. \(\rho ^0 + n\) events with additional activity in the forward detectors. In this sample the background fraction is found to be \(0.58\pm 0.07\). Fixing the DIFFVM normalisation by a fit to the ‘anti-LRG’ sample results in a background contribution of \(F_{bg}=0.29\pm 0.05\) in the main sample. Since the signal-to-background decomposition fit in this cross check gives a worse \(\chi ^2\), the nominal value \(F_{bg}=0.34\) is used for the cross section determinations. The difference to the \(F_{bg}\) value determined in the nominal analysis, as described above, is well covered by systematic uncertainty of the LRG condition efficiency.

3.6 Cross section determination and systematic uncertainties

The cross sections are measured for the kinematic ranges as defined in the rightmost column of Table 1. From the observed number of ep events, \(N_{data}\), the bin-integrated \(\gamma p\) cross section in bin i is calculated as

$$\begin{aligned} \sigma _{\gamma p}^i = \frac{1}{\Phi _{\gamma }} \frac{N_{data}^i - N_{bg}^i}{\mathcal{L}(A \cdot \epsilon )_i} \cdot C_{\rho }^i \end{aligned}$$
(11)

where \(N_{bg}^i\) is the expected diffractive dissociation background in bin i taking into account the overall normalisation fraction \(F_{bg}=0.34\), \(A \cdot \epsilon \) is the correction for detector acceptance and efficiency, \(\mathcal{L}\) is the integrated luminosity of the data, \(C_{\rho }\) is the extrapolation factor for the number of \(\rho ^0\) events from the \(M_{\pi \pi }\) measurement interval to the full \(\rho ^0\) mass range and \(\Phi _{\gamma }=0.1543\) is the value of the equivalent photon flux from Eq. (3) for the given \((W_{\gamma p},Q^2)\) range.Footnote 5 Since the statistics available does not allow for a reliable \(\rho ^0\) mass fit in every measurement bin, \(C_{\rho }^i\) is calculated using \(C_{\rho }=1.155\), obtained from the fit of the full sample and bin-dependent skewing correction factor derived from the fitted dependence of \(n_{RS}(p_{T,\rho }^2)\) in Eq. (9).

Several sources of experimental uncertainties are considered and their effects on the measured cross section are quantified. The systematic uncertainties on the cross section measurements are determined using MC simulations, by propagating the corresponding uncertainty through the full analysis chain. The individual systematic uncertainties are grouped into four categories below.

  • Detector related sources. The trigger efficiency is verified and tuned with the precision of \(3.4~\%\) using an independent monitoring sample. It is treated as correlated between different bins. The uncertainty due to the track finding and reconstruction efficiency in the central tracker is estimated to be \(1~\%\) per track [75] resulting in \(2~\%\) uncertainty in the cross section, taken to be correlated between bins. Several sources of uncertainties related to the measurement of the forward neutrons are considered. The uncertainty in the neutron detection efficiency which affects the measurement in a global way is \(2~\%\) [7]. The \(2~\%\) uncertainty on the absolute hadronic energy scale of the FNC [7] leads to a systematic error of \(1.1~\%\) for the \(x_L\)-integrated cross section and varying between 2 and \(19~\%\) in different \(x_L\) bins. 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 [7]. This results in an average uncertainty on the FNC acceptance determination of \(4.5~\%\) reaching up to \(10~\%\) for the \(p_{T,n}\) distribution. The systematics due to the exclusivity condition in the main part of the H1 detector is estimated to be \(2.1~\%\). It gets contributions from varying the LAr calorimeter noise cut between 400 MeV and 800 MeV (\(0.9~\%\)) and from the parameters of the algorithm connecting clusters with tracks (\(1.9~\%\)). This error influences only the overall normalisation. The uncertainty from the LRG condition is determined in the same manner as in the H1 inclusive diffraction analyses based on the large rapidity gap technique [76, 77]. It is further verified by comparing the cross sections obtained using different components of the forward detector apparatus for the LRG selection: FMD alone vs FMD\(+\)FTS vs FMD\(+\)Plug vs FMD\(+\)FTS\(+\)Plug. The resulting uncertainty is conservatively estimated to be \(9.0~\%\) affecting all bins in a correlated manner.

  • Backgrounds. Three different types of background are considered. Non-ep background is estimated from the shape of the z-vertex distribution and from the analysis of non-colliding proton bunches to be \((1.2 \pm 0.7)~\%\). Background originating from random coincidences between \(\rho ^0\) photoproduction events and neutrons from p-gas interactions amounts to \((1.0 \pm 0.2)~\%\). This results in \(2.2~\%\) background which was statistically subtracted in all distributions with an uncertainty of \(0.8~\%\). Non-\(\rho ^0\) background, as discussed in Sect. 3.4, has an uncertainty of \(0.7~\%\) and affects the overall normalisation only. Diffractive background to the DPP signal events (Sect. 3.5) is estimated with a precision of \(7.6~\%\). This is one of the largest individual uncertainties in the analysis. It is correlated between the bins.

  • MC model uncertainties. The uncertainty in the subtracted diffractive background due to the limited knowledge on \(\gamma p\) diffraction is evaluated by varying the \(M_Y\) and t dependencies in the DIFFVM simulation and the relative composition of diffractive channels within the limits allowed by previous HERA measurements. The resulting uncertainty is a part of the background subtraction systematics listed above. The systematic uncertainty of the MC correction factors for signal events is \(4.1~\%\), varying between \(1~\%\) and \(9~\%\) in different bins. It is evaluated from the difference between two versions of the POMPYT MC program with different \(p_T^2\) dependencies of the \(\rho ^0\) cross section, as described in Sect. 3.3. Here the uncertainty due to the POMPYT reweighting procedure is also accounted for.

  • Normalisation uncertainties. The uncertainty related to the \(\rho ^0\) mass fit, extrapolating from the measurement domain \(0.6 \le M_{\pi \pi } \le 1.1~\,\text{ GeV }\) to the full mass range \(2m_{\pi }<M_{\pi \pi }<M_{\rho }+5\Gamma _{\rho }\), which implies a correction factor of \( C_{\rho }=1.155\) on average in Eq. (11) with an uncertainty of \(1.6~\%\) due to fit errors. The integrated luminosity of the data sample is known with \(2.7~\%\) precision [46]. Together with other normalisation errors listed above the resulting total normalisation uncertainty amounts to \(4.4~\%\).

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 larger than the statistical uncertainties in most of the measurement bins.

The total systematic uncertainty for the integrated \(\gamma p\) cross section is \(14.6~\%\) including the global normalisation errors.

4 Results

Total, single- and double-differential photoproduction cross sections for the reaction \(\gamma p \rightarrow \rho ^0 n \pi ^+\) are measured in the kinematic range defined in Table 1. The photon-pion cross section, \(\sigma _{\gamma \pi }= \sigma (\gamma \pi ^+ \rightarrow \rho ^0\pi ^+)\), is extracted from the differential cross section \(\mathrm{d}\sigma _{\gamma p}/\mathrm{d}x_L\) using the pion flux [34, 35] integrated over the range \(p_{T,n}<0.2\) GeV. The results are summarised in Tables 2, 3, 4, 5, 6, 7, 8 and 9 and are shown in Figs. 6, 7, 8, 9, 10, 11, 12 and 13.

Fig. 6
figure 6

Differential cross section d\(\sigma _{\gamma p}/\mathrm{d}x_L\) in the range \(20<W_{\gamma p}<100\) GeV compared to the predictions based on different versions of the pion flux models. Top row Cross sections in the full FNC acceptance \(\theta _n<0.75\) mrad. Bottom row Cross sections for the OPE1 range, \(p_{T,n}<200\) MeV. Left-hand column Disfavoured versions of the pion fluxes, right-hand column pion fluxes compatible with the data. The data points are shown with statistical (inner error bars) and total (outer error bars) uncertainties, excluding an overall normalisation error of \(4.4~\%\). All predictions are normalised to the data

Table 2 Differential photoproduction cross sections \(\mathrm{d}\sigma _{\gamma p}/\mathrm{d}x_L\) for the exclusive process \(\gamma p \rightarrow \rho ^0 n \pi ^+\) in two regions of neutron transverse momentum and \(20\!<\!W_{\gamma p}\!<\!100\) GeV. The statistical, uncorrelated and correlated systematic uncertainties, \(\delta _{stat}\), \(\delta _{sys}^{unc}\) and \(\delta _{sys}^{cor}\) respectively, are given together with the total uncertainty \(\delta _{tot}\), which does not include the global normalisation error of \(4.4~\%\)
Table 3 Double differential photoproduction cross sections \(\mathrm{d^2}\sigma _{\gamma p}/\mathrm{d}x_L\mathrm{d}p_{T,n}^2\) in the range \(20\!<\!W_{\gamma p}\!<\!100\) GeV. The statistical, uncorrelated and correlated systematic uncertainties, \(\delta _{stat}\), \(\delta _{sys}^{unc}\) and \(\delta _{sys}^{cor}\) respectively, are given together with the total uncertainty \(\delta _{tot}\), which does not include the global normalisation error of \(4.4~\%\)
Table 4 The effective exponential slope, \(b_n\), obtained from the fit of double differential photoproduction cross sections \(\mathrm{d^2}\sigma _{\gamma p}/\mathrm{d}x_L\mathrm{d}p_{T,n}^2\) to a single exponential function in bins of \(x_L\). The first uncertainty represents the fit error from the statistical and uncorrelated systematic uncertainty and the second one is due to the correlated systematic uncertainty
Table 5 Energy dependence of the exclusive photoproduction of a \(\rho ^0\) meson associated with a leading neutron, \(\gamma p \rightarrow \rho ^0 n \pi ^+\). The first uncertainty is statistical and the second is systematic. The global normalisation uncertainty of \(4.4~\%\) is not included. \(\Phi _{\gamma }\) is the integral of the photon flux (3) in a given \(W_{\gamma p}\) bin
Table 6 Differential photoproduction cross section \(\mathrm{d}\sigma _{\gamma p}/\mathrm{d}\eta \) for the exclusive process \(\gamma p \rightarrow \rho ^0 n \pi ^+\) as a function of the \(\rho ^0\) pseudorapidity in the kinematic range \(0.35\!<\!x_L\!<\!0.95\), \(\theta _n\!<\!0.75\) mrad and \(20\!<\!W_{\gamma p}\!<\!100\) GeV. The statistical, uncorrelated and correlated systematic uncertainties, \(\delta _{stat}\), \(\delta _{sys}^{unc}\) and \(\delta _{sys}^{cor}\) respectively, are given together with the total uncertainty \(\delta _{tot}\), which does not include the global normalisation error of \(4.4~\%\)
Table 7 Differential photoproduction cross section \(\mathrm{d}\sigma _{\gamma p}/\mathrm{d}t^{\prime }\) for the exclusive process \(\gamma p \rightarrow \rho ^0 n \pi ^+\) as a function of the \(\rho ^0\) four-momentum transfer squared, \(t^{\prime }\), in the kinematic range \(0.35\!<\!x_L\!<\!0.95\), \(\theta _n\!<\!0.75\) mrad and \(20\!<\!W_{\gamma p}\!<\!100\) GeV. The statistical, uncorrelated and correlated systematic uncertainties, \(\delta _{stat}\), \(\delta _{sys}^{unc}\) and \(\delta _{sys}^{cor}\) respectively, are given together with the total uncertainty \(\delta _{tot}\), which does not include the global normalisation error of \(4.4~\%\)
Table 8 Exponential slopes, \(b_1\) and \(b_2\), and the ratio \(\sigma _1/\sigma _2\), obtained from the components of fit (14) to the differential cross section \(\mathrm{d}\sigma _{\gamma p}/\mathrm{d}t^{\prime }\) in bins of \(x_L\) and in bins of \(p_{T,n}^2\). The errors represent the statistical and systematic uncertainties added in quadrature
Table 9 Energy dependence of elastic \(\rho ^0\) photoproduction cross section on the pion, \(\gamma \pi ^+ \rightarrow \rho ^0 \pi ^+\), extracted in the one-pion-exchange approximation using OPE1 sample. The first uncertainty represents the full experimental error and the second is the model error coming from the pion flux uncertainty (see text). \(\Gamma _{\pi }(x_L)\) represents the value of the pion flux (5-6) integrated over the \(p_{T,n}<0.2\) GeV range, at a given \(x_L\)

4.1 \(\gamma p\) cross sections

The \(\gamma p\) cross section integrated in the domain \(0.35\!<\! x_L\!<\! 0.95\) and \(-t^{\prime }\!<\! 1\) GeV\(^2\) and averaged over the energy range \(20 \!<\! W_{\gamma p}\!<\! 100\) GeV is determined for two intervals of leading neutron transverse momentum:

$$\begin{aligned} \sigma (\gamma p \rightarrow \rho ^0 n \pi ^+ )= & {} (310 \pm 6_{\mathrm{stat}} \pm 45_\mathrm{sys})~ \mathrm{nb}\nonumber \\&\quad \hbox {for}\ p_{T,n}<x_L \cdot 0.69~\mathrm{GeV} \end{aligned}$$
(12)

and

$$\begin{aligned} \sigma (\gamma p \rightarrow \rho ^0 n \pi ^+ )= & {} (130 \pm 3_{\mathrm{stat}} \pm 19_\mathrm{sys})~ \mathrm{nb}\nonumber \\&\quad \mathrm{for}\ p_{T,n}<0.2~\mathrm{GeV}. \end{aligned}$$
(13)

Single differential cross sections as a function of \(x_L\) for these two regions are given in Table 2 and are shown in Fig. 6. The data are compared in shape to the predictions based on different models for the pion flux. Some models, like FMS [39] and NSSS [40] are disfavoured by the data and can be ruled out even on the basis of shape comparison alone. The other pion flux parametrisations: Bishari-0 [33], Holtmann [34, 35], KPP [37] and MST [38] are in good agreement with the data in both \(p_{T,n}\) ranges.

Additional constraints on the pion flux models could be provided by the dependence on t (or \(p_{T,n}^2\)) of the leading neutron. The double differential cross section \(\mathrm{d^2}\sigma _{\gamma p}/\mathrm{d}x_L\mathrm{d}p_{T,n}^2\) is measured, and the results are presented in Table 3 and Fig. 7. The bins are chosen such, that the data are not affected by the polar angle cut (see Fig. 2c). Although neither the t-, nor the \(p_{T,n}^2\)-dependence of the pion flux models are exactly exponential they can be approximated by a simple exponent in many cases. Such an approximation has been used already in other analyses [2, 6]. The \(p_{T,n}^2\)-distributions measured here for fixed \(x_L\) are compatible with an exponential shape within the statistical and uncorrelated systematic errors. Therefore, the same approach is used here. The cross sections are fitted by a single exponential function \(e^{-b_n(x_L)p_{T,n}^2}\) in each \(x_L\) bin. The quality of the fits is good, with \(P(\chi ^2)=0.35 \div 0.60\). The results are presented in Table 4 and Fig. 8. The measured b-slopes are compared to those obtained from several pion flux parametrisations. Despite of the large experimental uncertainties none of the models is able to reproduce the data.Footnote 6 A possible reason for this discrepancy could be the effect of energy-momentum conservation affecting the proton vertex in this exclusive reaction more strongly than in inclusive production of a leading neutron in which an apparent factorisation of the proton vertex has been observed. Another explanation [21, 89] could be absorptive corrections which modify the t dependence of the amplitude, leading to an increase of the effective b-slope at large \(x_L\) as compared to the pure OPE model without absorption.

Fig. 7
figure 7

Double differential cross section \(\mathrm{d^2}\sigma _{\gamma p}/\mathrm{d}x_L\mathrm{d}p_{T,n}^2\) of neutrons in the range \(20<W_{\gamma p}<100\) GeV fitted with single exponential functions. The cross sections in different \(x_L\) bins j are scaled by the factor \(3^j\) for better visibility. The binning scheme is shown in Fig. 2c. The data points are shown with statistical (inner error bars) and total (outer error bars) uncertainties excluding an overall normalisation error of \(4.4~\%\)

Fig. 8
figure 8

The exponential slopes fitted through the \(p_T^2\) dependence of the leading neutrons as a function of \(x_L\). The inner error bars represent statistical errors and the outer error bars are statistical and systematic errors added in quadrature. The data points are compared to the expectations of several parametrisations of the pion flux within the OPE model

Fig. 9
figure 9

Cross section of the reaction \(\gamma p \rightarrow \rho ^0 n \pi ^+\) as a function of \(W_{\gamma p}\) compared to the prediction from POMPYT MC program, which is normalised to the data. The dashed curve represents the Regge motivated fit \(\sigma \propto W^{\delta }\) with \(\delta = -0.26 \pm 0.06_{\mathrm{stat}}\pm 0.07_\mathrm{sys}\). The data points are shown with statistical (inner error bars) and total uncertainties (outer error bars) excluding an overall normalisation error of \(4.4~\%\)

The energy dependence of the reaction \(\gamma p \rightarrow \rho ^0 n \pi ^+\) is presented in Table 5 and in Fig. 9. The cross section drops with \(W_{\gamma p}\) in contrast to the POMPYT MC expectation, where the energy dependence is driven by Pomeron exchange alone. A Regge motivated power law fit to the data, \(\sigma _{\gamma p}(W_{\gamma p}) \propto W_{\gamma p}^{\delta }\), yields \(\delta = -0.26 \pm 0.06_{stat} \pm 0.07_{sys}\). The difference in the energy dependence in data and MC is also reflected in the pseudorapidity distribution of the \(\rho ^0\) meson, which is given in Table 6 and shown in Fig. 10.

Fig. 10
figure 10

The differential cross section d\(\sigma _{\gamma p}/\mathrm{d}\eta \) as a function of pseudorapidity of the \(\rho ^0\) meson compared to the prediction from the POMPYT MC program, which is normalised to the data. The inner error bars represent statistical errors and the outer error bars are total errors excluding an overall normalisation uncertainty of \(4.4~\%\)

Fig. 11
figure 11

Differential cross section d\(\sigma _{\gamma p}/\mathrm{d}t^{\prime }\) of \(\rho ^0\) mesons fitted with the sum of two exponential functions. The inner error bars represent statistical and uncorrelated systematic uncertainties added in quadrature and the outer error bars are total uncertainties, excluding an overall normalisation error of \(4.4~\%\)

Fig. 12
figure 12

The two exponential slopes, \(b_1\) and \(b_2\), obtained by fitting the \(t^{\prime }\) dependence of the \(\rho ^0\) mesons (left) and the relative contribution of the two exponents to the overall cross section of the reaction \(\gamma p \rightarrow \rho ^0 n \pi ^+\) (right) as a function of \(x_L\) (top) and \(p_{T,n}^2\) (bottom). The error bars represent statistical and systematic uncertainties added in quadrature. Horizontal lines with error bands show the corresponding average values for the full ranges of \(0.35<x_L<0.95\) and \(0<p_{T,n}^2< 0.3~\mathrm {GeV}^2 \)

Fig. 13
figure 13

Elastic cross section, \(\sigma _{\gamma \pi }^\mathrm{el} \equiv \sigma ({\gamma \pi ^+} \rightarrow \rho ^0\pi ^+)\), extracted in the one-pion-exchange approximation as a function of the photon-pion energy, \(W_{\gamma \pi }\). The inner error bars represent the total experimental uncertainty and the outer error bars are experimental and model uncertainties added in quadrature, where the model error is due to pion flux uncertainties. The dark shaded band represents the average value for the full \(W_{\gamma \pi }\) range as given in Eq. (16)

Finally, the cross section as a function of the four-momentum transfer squared of the \(\rho ^0\) meson, \(t^{\prime }\), is given in Table 7 and presented in Fig. 11. It exhibits the very pronounced feature of a strongly changing slope between the low-\(t^{\prime }\) and the high-\(t^{\prime }\) regions. The fit is performed to the sum of two exponential functions:

$$\begin{aligned} \frac{\mathrm{d}\sigma _{\gamma p}}{\mathrm{d}t^{\prime }} = a_1 e^{b_1 t^{\prime }} + a_2 e^{b_2 t^{\prime }} \end{aligned}$$
(14)

and yields the following slope parameters:

$$\begin{aligned}&b_1 = (25.72 \pm 3.22_{unc} \pm 0.26_{cor})~\mathrm{GeV}^{-2}; \nonumber \\&b_2 = (3.62 \pm 0.30_{unc} \pm 0.10_{cor})~\mathrm{GeV}^{-2} \end{aligned}$$
(15)

where the first errors include statistical and uncorrelated systematic uncertainties and the second errors are due to correlated systematic uncertainties. In a geometric picture, the large value of \(b_1\) suggests that for a significant part of the data \(\rho ^0\) mesons are produced at large impact parameter values of order \(\langle r^2\rangle = 2b_1\!\cdot \!(\hbar c)^2 \simeq 2~\mathrm{fm}^2 \approx (1.6 R_\mathrm{p})^2\). In other words, photons find pions in a cloud which extends far beyond the proton radius. The small value of \(b_2\) corresponds to a target size of \({\sim } 0.5\) fm. In the DPP interpretation [21, 24, 25] the observed behaviour is a consequence of the interference between the amplitudes corresponding to the diagrams a, b and c in Fig. 1, leading to a slope dependence on the invariant mass of the \((n\pi ^+)\) system produced at the proton vertex. Since the forward pion is not detected in this analysis the \((n\pi ^+)\) invariant mass cannot be determined with sufficient precision, which prevents explicit measurement of the b(m) dependence.

In order to investigate the presence of a possible factorisation between the proton and the photon vertices, the \(t^{\prime }\) distribution is studied in bins of \(x_L\). The result of the fit by Eq. (14) with \(x_L\) dependent parameters \(a_i(x_L), b_i(x_L)\) is presented in Table 8 and in Fig. 12 in comparison with the values given in Eq. (15) for the full \(x_L\) range. Also the evolution with \(x_L\) of the ratio of two components, \(\sigma _1/\sigma _2\), where \(\sigma _i = \frac{a_i}{b_i}(1-e^{-b_i})\), is shown. Given the large experimental uncertainties no strong conclusion about factorisation of the two vertices can be drawn.

4.2 \(\gamma \pi \) cross section

The pion flux models compatible with the data in shape of the \(x_L\) distribution are used to extract the photon-pion cross sections from \(\mathrm{d}\sigma /\mathrm{d}x_L\) in the OPE approximation. The results are presented in Table 9 and in Fig. 13. As a central value the Holtmann flux [34, 35] is used, and the largest difference to the other three predictions [33, 37, 38] provides an estimate of the model uncertainty which is \({\sim }19~\%\) on average. From the total \(\gamma p\) cross section in Eq. (13) and using the pion flux (5) and (6) integrated in \(x_L\) and \(p_{T,n}\), \(\Gamma _{\pi } = 0.056\), the cross section of elastic photoproduction of \(\rho ^0\) on a pion target is determined at an average energy \(\langle W_{\gamma \pi }\rangle \simeq 24\) GeV:

$$\begin{aligned} \sigma (\gamma \pi ^+\rightarrow \rho ^0\pi ^+) = (2.33\pm 0.34 (\mathrm exp) ^{+0.47}_{-0.40} (\mathrm model))~\mu \mathrm b,\nonumber \\ \end{aligned}$$
(16)

where the experimental uncertainty includes statistical, systematic and normalisation errors added in quadrature, while the model error is due to the uncertainty in the pion flux integral obtained for the different flux parametrisations compatible with our data.

Theoretical studies of leading neutron production in ep collisions [34, 35, 37] suggest that in addition to the pion exchange process other processesFootnote 7 may contribute at \(10{-}20~\%\) level. To suppress these contributions it is recommended to perform cross section measurements in the ‘OPE safe’ phase space region: low \(p_{T,n}\) and high \(x_L\). In order to investigate a possible influence of non-OPE contributions the extraction of the photon-pion cross section is repeated for two additional regions, in which the validity of pure OPE is assumed. The cross sections for the full FNC acceptance range \((\theta _n<0.75~\mathrm{mrad},~0.35<x_L<0.95)\) and for the OPE2 sub-sample \((p_{T,n}<200~\mathrm{MeV},~0.65<x_L<0.95)\) together with the value (16) obtained for the OPE1 sample are presented in Table 10. The values of \(\sigma (\gamma \pi ^+\rightarrow \rho ^0\pi ^+)\) extracted in these three different phase space regions agree well within the experimental errors. Thus no evidence for an extra contribution beyond the OPE is found in the full FNC acceptance range for the exclusive reaction studied here.

Table 10 Cross section of elastic \(\rho ^0\) photoproduction on the pion, \(\gamma \pi ^+ \rightarrow \rho ^0 \pi ^+\), extracted in the one-pion-exchange approximation using three different samples: full sample, OPE1 and OPE2. The first uncertainty represents the full experimental error and the second is the model error coming from the pion flux uncertainty (see text). \(\Gamma _{\pi }\) represents the value of the pion flux (5) and (6) integrated over the corresponding \((x_L,p_{T,n})\) range

Taking a value of \(\sigma (\gamma p \rightarrow \rho ^0 p) = (9.5 \pm 0.5)~\mu \mathrm{b}\) at the corresponding energy \(\langle W \rangle =24\) GeV, which is an interpolation between fixed target and HERA measurements (see e.g. figure 10 in [56]), one obtains for the ratio \(r_\mathrm{el} = \sigma _\mathrm{el}^{\gamma \pi }/\sigma _\mathrm{el}^{\gamma p} = 0.25 \pm 0.06\). A similar ratio, but for the total cross sections at \(\langle W \rangle =107\) GeV, has been estimated by the ZEUS collaboration as \(r_\mathrm{tot} = \sigma _\mathrm{tot}^{\gamma \pi }/\sigma _\mathrm{tot}^{\gamma p} = 0.32 \pm 0.03\) [2]. Both ratios are significantly smaller than their respective expectations, based on simple considerations. For \(r_\mathrm{tot}\), a value of 2 / 3 is predicted by the additive quark model [7880], while \(r_\mathrm{el} = (\frac{b_{\gamma p}}{b_{\gamma \pi }}) \cdot (\sigma _\mathrm{tot}^{\gamma \pi }/\sigma _\mathrm{tot}^{\gamma p})^2 = 0.57 \pm 0.03\) can be deduced by combining the optical theorem, the eikonal approach [81] relating cross sections with elastic slope parameters [82] and the data on \(pp, \pi ^+p\) [8385] and \(\gamma p\) [57] elastic scattering. Such a suppression of the cross section is usually attributed to rescattering, or absorptive corrections [8689], which are essential for leading neutron production. For the exclusive reaction \(\gamma p \rightarrow \rho ^0 n \pi ^+\) studied here this would imply an absorption factor of \(K_{abs} = 0.44 \pm 0.11\). It is interesting to note, that this value is similar to the somewhat different, but conceptually related damping factor in diffractive dijet photoproduction, the rapidity gap survival probability, \(\langle S^2 \rangle \simeq 0.5\), which has been determined by the H1 collaboration [9092].

5 Summary

The photoproduction cross section for exclusive \(\rho ^0\) production associated with a leading neutron is measured for the first time at HERA. The integrated \(\gamma p\) cross section in the kinematic range \(20<W_{\gamma p}<100\) GeV, \(0.35<x_L<0.95\) and \(\theta _n<0.75\) mrad is determined with \(2~\%\) statistical and \(14.6~\%\) systematic precision. The elastic photon-pion cross section, \(\sigma (\gamma \pi ^+ \rightarrow \rho ^0\pi ^+)\), at \(\langle W_{\gamma \pi }\rangle = 24\) GeV is extracted in the one-pion-exchange approximation.

Single and double differential \(\gamma p\) cross sections are measured. The differential cross section d\(\sigma /\mathrm{d}t^{\prime }\) shows a behaviour typical for exclusive double peripheral exchange processes.

The differential cross sections for the leading neutron are sensitive to the pion flux models. While the shape of the \(x_L\) distribution is well reproduced by most of the pion flux parametrisations, the \(x_L\) dependence of the \(p_T\) slope of the leading neutron is not described by any of the existing models. This may indicate that the proton vertex factorisation hypothesis does not hold in exclusive photoproduction, e.g. due to large absorptive effects which are expected to play an essential rôle in soft peripheral processes. The estimated cross section ratio for the elastic photoproduction of \(\rho ^0\) mesons on the pion and on the proton, \(r_\mathrm{el} = \sigma _\mathrm{el}^{\gamma \pi }/\sigma _\mathrm{el}^{\gamma p} = 0.25 \pm 0.06\), suggests large absorption corrections, of the order of \(60~\%\), suppressing the rate of the studied reaction \(\gamma p \rightarrow \rho ^0 n \pi ^+\).