Gluon saturation and Feynman scaling in leading neutron production

In this paper we extend the color dipole formalism to the study of leading neutron production in $e + p \rightarrow e + n + X$ collisions at high energies and estimate the related observables, which were measured at HERA and may be analysed in future electron-proton ($ep$) colliders. In particular, we calculate the Feynman $x_F$ distribution of leading neutrons, which is expressed in terms of the pion flux and the photon-pion total cross section. In the color dipole formalism, the photon-pion cross section is described in terms of the dipole-pion scattering amplitude, which contains information about the QCD dynamics at high energies and gluon saturation effects. We consider different models for the scattering amplitude, which have been used to describe the inclusive and diffractive $ep$ HERA data. Moreover, the model dependence of our predictions with the description of the pion flux is analysed in detail. We show that the recently released H1 leading neutron spectra can be reproduced using the color dipole formalism and that these spectra could help us to observe more clearly gluon saturation effects in future $ep$ colliders.


I. INTRODUCTION
In high energy collisions the outgoing baryons which have large fractional longitudinal momentum (x L ≥ 0.3) and the same valence quarks (or at least one of them) as the incoming particles are called leading particles (LP). The momentum spectra of leading particles have been measured already some time ago [1]. Very recently, high precision data on leading neutrons produced in electron-proton reactions at HERA at high energies became available [2].
Leading neutron production is a very interesting process. In spite of more than ten years of intense experimental and theoretical efforts, the x L (Feynman momentum) distribution of leading neutrons remains without a satisfactory theoretical description [3][4][5][6][7][8][9][10][11][12]. Monte Carlo studies, using standard deep inelastic scattering (DIS) generators show [2,13] that these processes have a rate of neutron production a factor of three lower than the data and produce a neutron energy spectrum with the wrong shape, peaking at values of x L below 0.3. In order to fit the data the existing models need to combine different ingredients including pion exchange, reggeon exchange, baryon resonance excitation and decay and independent fragmentation [8][9][10][11][12]. Moreover, the incoming photon interacts with the pion emitted by proton and then rescatters, interacting also with the emerging neutron and giving origin to significant absorptive corrections, which are difficult to calculate [11,12]. As it can be seen in Fig. 1, this process is essentially composed by a soft pion (or reggeon) emission and by the subsequent photon-pion interaction at high energies. Pion emission has been studied for a long time and according to the traditional wisdom it can be described by a simple interaction Lagrangian of the form gψγ 5 πψ, where ψ and π are the nucleon and pion field respectively. The corresponding pionnucleon amplitude must be supplemented with a form factor, which represents the extended nature of hadrons and at the same time regularizes divergent integrals. The precise functional form of the form factor is very hard (if not impossible) to extract from first principle calculations. We must then resort to models. Very recently [14,15] it has been argued that the Lagrangian which is more consistent with the chiral symmetry requirements is the one with a pseudovector couplings. One of the goals of the present study is to test in a new context the better founded splitting function, f (y), derived in [14]. Additionally, forward hadron production is very important also for cosmic ray physics, where highly energetic protons reach the top of the atmosphere and undergo successive high energy scatterings off light nuclei in the air. In each of these collisions, a projectile proton (the leading baryon) looses energy, creating showers of particles, and goes to the next scattering. The interpretation of cosmic data depends on the accurate knowledge of the leading baryon momentum spectrum and its energy dependence. A crucial question of practical importance is the existence or non-existence of the Feynman scaling, which says that the x L spectra of secondaries are energy independent. In cosmic ray applications we are sensitive essentially to the large x L region (the fragmentation region), which probes the low Bjorken-x component of the target wave function. In this kinematical range nonlinear effects are expected to be present in the description of the QCD dynamics (for recent reviews see [16]), associated to the high parton density. The state-of-art framework to treat QCD at high energies is the Color Glass Condensate (CGC) formalism [17], which predicts gluon saturation at small-x, with the evolution with the energy being described by an infinite hierarchy of coupled equations for the correlators of Wilson lines -the Balitsky-JIMWLK hierarchy [16]. In the mean field approximation, this set of equations can be approximated by the Balitsky-Kovchegov (BK) equation [18,19], which describes the evolution of the dipole-target scattering amplitude with the rapidity Y = ln(1/x). In this paper we propose to treat the leading neutron production at HERA using the color dipole formalism, which is able to describe the inclusive and diffractive HERA ep data. Our goal is to extend this successful formalism, which have its main parameters well determined, for leading neutron production. As a consequence, our predictions are free parameter, depending only from the choices for the models of the pion flux and dipole scattering amplitude. Moreover, the use of the color dipole formalism allows to estimate the contribution of the gluon saturation effects for the leading neutron production in the kinematical range which was probed by HERA and which will probed in future electron-proton colliders. Finally, it allows to investigate the relation between the Feynman scaling (or its violation) and the description of the QCD dynamics at high energies. It is important to emphasize that, in the near future, Feynman scaling will be investigated experimentally at the LHC by the LHCf Collaboration [20][21][22]. This paper is organized as follows. In the next Section we present a brief review of the leading neutron production in ep collisions as well as the different models for the pion flux are discussed. Moreover, the treatment of the process using the color dipole formalism is presented and the main assumptions are analysed. In Section III we analyse the dependence of our predictions in the pion flux and in the scattering amplitude. A comparison with the recent H1 data is presented and the Feynman scaling is analysed. Finally, in Section IV we summarize our main conclusions.

A. The cross section
Let us review the main formulas of leading neutron production. At high energies, this scattering can be seen as a set of three factorizable subprocesses (See Fig. 1): i) the photon fluctuates into a quark-antiquark pair (the color dipole), ii) the color dipole interacts with the pion, present in the wave function of the incident proton, and iii) the leading neutron is formed. The differential cross section reads: where Q 2 is the virtuality of the exchanged photon,Ŵ is the center-of-mass energy of the virtual photon-pion system. It can be written asŴ 2 = (1 − x L ) W 2 , where W is the center-of-mass energy of the virtual photon-proton system. As it can be seen in Fig. 1, x L is the proton momentum fraction carried by the neutron and t is the square of the four-momentum of the exchanged pion. In terms of the measured quantities x L and transverse momentum p T , the pion virtuality is: The flux of virtual pions emitted by the proton is represented by f π/p and σ γ * π (Ŵ 2 , Q 2 ) is the cross section for the interaction between the virtual-photon and the virtual-pion at center-of-mass energyŴ .

B. The pion flux
The pion flux f π/p (x L , t) (also called sometimes pion splitting function) is the virtual pion momentum distribution in a physical nucleon (the bare nucleon plus the "pion cloud"). It was first calculated in the early studies [23] of deep inelastic scattering (DIS), where a pseudoscalar nucleon-pion-nucleon vertex was added to the standard DIS diagram. These early calculations were further refined in Refs. [24] and even extended to the strange and charm sector [25]. Since pion emission is a nonperturbative process, once we depart from the single nucleon state, we should consider a whole tower of meson-baryon states, having to deal with a series for which there is no rigorous truncation scheme. This led some authors to use light cone models [26] for the pion momentum distribution, where the dynamical origin of the pion is not mentioned and one tries to determine phenomenologically the relative weight of the higher Fock states.
In all the calculations of the pion flux a form factor was introduced to represent the non pointlike nature of hadrons and hadronic vertices, which contain a cut-off parameter determined by fitting data. The most frequently used parametrizations of the pion flux [3][4][5][6][7][8][9][10][11][12] have the following general form: where g 2 pπp /(4π) = 14.4 is the π 0 pp coupling constant, m π is the pion mass and α(t) will be defined below. The form factor F (x L , t) accounts for the finite size of the nucleon and pion. We will consider the following parametrizations of the form factor: , α(t) = 0 (4) from Ref. [4], where R = 0.6 GeV −1 .
from Ref. [3], where α π (t) ≃ t (with t in GeV 2 ) is the Regge trajectory of the pion.
also from Ref. [6], where Λ d = 1.2 GeV. In the case of the more familiar exponential (4), monopole (7) and dipole (8) forms factors, the cut-off parameters have been determined by fitting low energy data on nucleon and nuclear reactions and also data on deep inelastic scattering and structure functions.
In all analyses of pion cloud effects in DIS the authors have used pseudoscalar (PS) coupling, which is, by itself, inconsistent with chiral symmetry. More recently, in [14] the authors used an effective chiral Lagrangian for the interaction of pions and nucleons consistent with chiral symmetry. Unlike the earlier chiral effective theory calculations which only computed the light-cone distributions of pions or considered the non analytic behaviour of their moments to lowest order in the pion mass, in [14] the authors have computed the complete set of diagrams relevant for DIS from nucleons dressed by pions resulting from their Lagrangians, without taking the heavy baryon limit. In particular, they have demonstrated explicitly the consistency of the computed distribution functions with electromagnetic gauge invariance. In [15] they applied the previously computed pion distribution to the study of thed −ū asymmetry in the nucleon. For our purposes, the relevant contribution of the pion momentum distribution reads [14]: where y = 1 − x L , g A = 1.267 is the nucleon axial charge, f π = 93 MeV and Λ c = 0.2 GeV. In order to obtain the final cross section we multiply (9) by 2 (see [14,15]) and insert it into (1) after integrating the latter over t (or over In what follows, we shall use the pion fluxes listed above, in equations (4)- (9), denoting them by f 1 , f 2 , ... f 6 , respectively. The main purpose of our calculation will be to show that the dipole approach gives a good description of data and the pion flux is just one element of the calculation. Therefore we will make no effort to choose one particular form. Nevertheless, we note that one important phenomenological constraint that these pion fluxes must satisfy is to reproduce thed −ū asymmetry in the proton sea measured by the E866 Collaboration [27]. Among the pion fluxes mentioned above only (8) (see [28]) and (9) have been confronted with these data. As it will be seen later, the results discussed here are very sensitive to the choice of the pion flux and hence leading neutron spectra may be used to constrain its shape.
Before closing this subsection we would like to mention that the diagram shown in Fig.1 represents only the dominant contribution to leading neutron production. Other isovector meson exchanges, such as ρ or a 2 , can also contribute to the leading neutron spectrum. Moreover the p → ∆ transition can also contribute to neutron production through the subsequent decay ∆ → nπ. Theoretical studies show that processes other than direct pion exchange are expected to give a contribution of less than 25 % of the cross section [4][5][6][7]29]. The effect of these backgrounds to the one pion exchange is to increase the rate of neutron production. However this effect is partially compensated by the absorptive rescattering of the neutron, which decreases the neutron rate by approximately the same amount [11,12].

C. The photon-pion cross section
In order to obtain the photon-pion cross section we will use the color dipole formalism, as usually done in high energy deep inelastic scattering off a nucleon target. In this formalism, the cross section is factorized in terms of the photon wave functions Ψ, which describes the photon splitting in a qq pair, and the dipole-pion cross section σ dπ . We have is the scaled Bjorken variable and the variable r defines the relative transverse separation of the pair (dipole). Moreover, we have that the photon wave functions are given by for a longitudinally (L) and transversely (T) polarized photon, respectively. In the above expressions ǫ 2 = z(1 − z)Q 2 + m 2 f , K 0 and K 1 are modified Bessel functions and the sum is over quarks of flavour f with a corresponding quark mass m f . As usual z stands for the longitudinal photon momentum fraction carried by the quark and 1 − z is the longitudinal photon momentum fraction of the antiquark. The main input in the calculations of σ γ * π is the dipole-pion cross section. In what follows, for simplicity, we will assume the validity of the additive quark model, which allows us to relate σ dπ with the dipole-proton cross section, usually probed in the typical inclusive and exclusive processes at HERA. Basically, we will assume that σ dπ (x, r) = 2 3 · σ dp (x, r) This is assumption is supported by the study of the pion structure function in the low x regime presented in [30]. It also gives a good description of the previous ZEUS leading neutron spectra, as shown in [9,10]. On the other hand, the direct application of (1) to HERA photoproduction data [31] leads to the result σ tot γπ /σ tot γp = 0.32 ± 0.03, which is factor 2 lower than the ratio given above. This subject certainly deserves more detailed studies in the future. For our purposes, the use of relation (14) allows us to estimate σ dπ without additional free parameters. In the eikonal approximation the dipole-proton cross section σ dp is given by: where N p (x, r, b) is the imaginary part of the forward amplitude for the scattering between a small dipole (a colorless quark-antiquark pair) and a dense hadron target, at a given rapidity interval Y = ln(1/x). The dipole has transverse size given by the vector r = x − y, where x and y are the transverse vectors for the quark and antiquark, respectively, and impact parameter b = (x + y)/2. As mentioned in the introduction, at very high energies (and very low x) the evolution with the rapidity Y of N p (x, r, b) is given by the Balitsky-Kovchegov (BK) equation [18,19] assuming the translational invariance approximation, which implies N p (x, r, b) = N p (x, r)S(b) and σ dp (x, r) = σ 0 · N p (x, r), with the normalization of the dipole cross section (σ 0 ) being fitted to data. Alternatively, we can describe the scattering amplitude N p (x, r) using phenomenological models based on saturation physics constructed taking into account the analytical solutions of the BK equation which are known in the low and high density regimes. The main advantage in the use of phenomenological models is that we can easily to compare the linear and nonlinear predictions, which is useful to determine the contribution of the saturation effects for the process under analysis. In what follows we will consider as input the phenomenological models proposed in Refs. [32][33][34]. In particular, the IIMS model, proposed in Ref. [33] and updated in [34], was constructed so as to reproduce two limits of the LO BK equation analytically under control: the solution of the BFKL equation for small dipole sizes, r ≪ 1/Q s (x), and the Levin-Tuchin law for larger ones, r ≫ 1/Q s (x). In the updated version of this parametrization [34], the free parameters were obtained by fitting the new H1 and ZEUS data. In this parametrization the forward dipole-proton scattering amplitude is given by where a and b are determined by continuity conditions at rQ s (x) = 2, γ s = 0.7376, κ = 9.9, N 0 = 0.7 and Q s is the saturation scale given by: with x 0 = 1.632 × 10 −5 , λ = 0.2197, Q 2 0 = 1.0 GeV 2 . The first line of Eq. (16) describes the linear regime whereas the second one includes saturation effects. On the other hand, in the GBW model [32], the dipole-proton scattering amplitude is given by with Q 2 0 = 1.0 GeV 2 , x 0 = 3 × 10 −4 and λ = 0.288. Finally, we will also use the dipole-proton cross section estimated with the help of the DGLAP analysis of the gluon distribution, which is given by [35]: where xg(x, Q 2 ) is the target gluon distribution, for which we use the CTEQ6 parametrization [36]. The above expression represents the linear regime of QCD and is a baseline for comparison with the nonlinear predictions.

III. RESULTS
As it was seen in the previous sections, the leading neutron spectrum depends essentially on two main ingredients: the pion flux and the dipole-pion cross section. The latter is very sensitive to the value of the involved Bjorken x, or, in our case,x. In Fig. 2 we showx as a function of x L . From the figure we can conclude that at the highest values of the present photon-hadron energies (W ) and at the lowest values of Q 2 we enter deeply in the low x domain. This will be even more so if measurements can be carried out at higher energies, but already at the available energies we can see that the leading neutron spectrum receives a contribution from the kinematical range where the nonlinear effects are expected to be present and hence it is interesting to investigate their influence.
We address now the dependence of our predictions for the leading neutron spectra on the models of the pion flux and of the dipole-pion scattering amplitude. It is important to emphasize that the spectrum is proportional to the product of these two quantities. Consequently, in what follows we will initially assume a given model for one of these quantities and analyse the dependence on the other quantity. Having done that, we will choose one combination, calculate the normalized cross section (1/σ DIS )dσ/dx L and compare with the experimental data. The normalization was taken from Ref. [37] and is given by: where c = 0.18, β = d ln(Q 2 /Λ 2 0 ) with d = 0.0481 and Λ 0 = 0.292 GeV. In Fig. 3 (a) we consider the IIMS model for the scattering amplitude and estimate the leading neutron spectrum using the different models of the pion flux discussed in the previous section. We observe that the behavior at medium and small values of x L is strongly dependent on the choice of the model. Using the cut-off Λ c in the range defined in Ref. [15], our results suggest that model f 6 is disfavored. In Fig. 3 (b) we consider the f 2 model for the pion flux and estimate the spectra using different models of the scattering amplitude. We can see that the predictions have similar behavior in the transition between the small and large-x L regimes, but with the magnitude being dependent on the description of the QCD dynamics. In the previous analysis we have considered Q 2 = 53 GeV 2 . However, the experimental data were obtained for a wide range of virtualities, which implies that there is some arbitrariness present in the choice of the value of Q 2 used in our calculations. In order to estimate the uncertainty present in this choice, in Fig. 3 (c) we consider the f 2 and IIMS models and calculate the spectrum for different values of the photon virtuality. We observe that our predictions are more compatible with the experimental data if larger values of Q 2 are assumed.
So far we have assumed that the photon hits only the pion, in a type of "impulse aproximation". However it has been shown in [12] that very often the photon hits also the neutron, specially in the low Q 2 domain, where the photon has low resolving power. In these cases the extra interactions generate the so called absorptive corrections, which can be estimated with models. In the case of leading neutron production, it has been shown [12] that the corrections are not very large and affect almost uniformly all the x L spectrum. Using the results of [12], we shall here simply multiply our spectra by the absorptive correction factor, K abs = 0.7.
In Fig. 4 we give our description of data, which is quite reasonable. All the parameters contained in the dipole cross section have been already fixed by the analysis of other DIS data from HERA. The pion flux has been fixed by fitting data of low energy hadronic collisions. As metioned above, it has been argued that the factor 2 3 predicted by the additive quark model to relate the photon-pion and photon-proton cross sections could be smaller, closer to 1 3 . Using this different factor, we would then underpredict the data in the low x L region. In principle, this would be acceptable since in this region other processes are expected to play a role. As already emphasized, this subject deserves more detailed studies.
Finally, we analyse the Feynman scaling in the leading neutron spectra and the contribution of nonlinear effects for this process. As mentioned before this process is surprisingly sensitive to low x physics. In particular, the large x L (x L > 0.3) we observe the transition from the large to small x domain. It is therefore interesting to check whether gluon saturation effects are already playing a significant role. As it can be seen from Eqs. (1) and (10), all the energy dependence is contained in the dipole forward scattering amplitude N . Therefore, we can estimate the contribution of the nonlinear effects comparing the results obtained with saturation model, e.g. the IIMS model given in (16), with those obtained with the linear model, Eq. (19). Moreover, the theoretical expectation can be obtained using the GBW model for the scattering amplitude, Eq. (18). In the linear limit, when the dipole radius is very small (or equivalently Q 2 is very large) or the saturation scale is very small (and hence the energy is not very high), we can expand the exponent and obtain Consequently, in this regime we see that the leading neutron x L spectrum will depend on W . In a complementary way, in the nonlinear limit, when the dipole radius is very large (or equivalently Q 2 is very small) or the saturation scale is very large (and hence the energy is very high), we obtain which is energy independent. One could argue that the mere inspection of Eq. (1) would suggest that at some asymptotically high energy the photon-pion cross section would reach some "black disk" limit and the energy dependence would disappear. We would like to emphasize that the information contained in (21) and (22) is much richer and indicates the route through which the asymptotic limit is reached and the role played by nonlinear effects. These expectations can be compared with those obtained using the IIMS and DGLAP models for the dipole-pion cross section. In Fig. 5 (a) we show the spectra obtained in a purely linear approach. As expected we see a noticeable energy dependence. In contrast, the nonlinear predictions presented in Fig. 5 (b) show a remarkable suppression of the energy dependence at low values of Q 2 , consistent with the expectations. These results indicate that the Feynman scaling (and how it is violated) can be directly related to the QCD dynamics at small-x.

IV. SUMMARY
In this work we have studied leading neutron production in e + p → e + n + X collisions at high energies and have calculated the Feynman x L distribution of these neutrons. The differential cross section was written in terms of the pion flux and the photon-pion total cross section. We have proposed to describe this process using the color dipole formalism and, assuming the validity of the additive quark model, we have related the dipole-pion with the well determined dipole-proton cross section. In this formalism we have been able to estimate the dependence of the predictions on the description of the QCD dynamics at high energies as well as the contribution of the gluon saturation effects for the leading neutron production. With the parameters strongly constrained by other phenomenological information, we were able to reproduce the recently released H1 leading neutron spectra. One of our most interesting conclusions is that leading neutron spectra can be used to probe the low x content of the pion target and hence it is a new observable where we can look for gluon saturation effects. At higher energies this statement will become more valid. Moreover saturation physics provides a precise route to Feynman scaling, which will eventually happen at higher energies. These conclusions motivate more detailed studies, in particular in the subjects which are the main source of theoretical uncertainties in our calculations: i) the validity of the additive quark model; ii) the factor K abs accounting for the absorptive corrections, which are model dependent; iii) the precise form of the pion flux; iv) the precise form of the dipole cross section. It is important to emphasize that none of these numbers or functions is free. On the contrary, they are subject to severe constraints from other experimental information and the freedom to choose them will be further reduced in the future.