QCD analysis of forward neutron production in DIS

We consider forward neutron production in DIS within fracture functions formalism. By performing a QCD analysis of available data we extract proton-to-neutron fracture functions exploiting a method which is in close relation with the factorisation theorem for this class of processes.


I. INTRODUCTION
In hadronic collisions a portion of the produced particle spectrum is characterised by hadrons carrying a sizeable fraction of the available energy and produced at small polar angle with respect to the collision axis. It is phenomenologically observed that for such hadrons their valence-parton composition is almost or totally conserved with respect to the one of initalstate hadrons [1]. Such semi-inclusive processes are instrumental to study the scaling hypothesis of forward hadron production cross sections [2,3] and give insight on non-perturbative aspects of QCD dynamics in high energy collisions. These hadrons, in fact, are produced at very small transverse momenta with respect to the collision axis, a regime where perturbative techniques can not be applied. Quite interestingly, forward particle production has also been observed in processes which involve point-like probes in lepton-hadron interactions, such as Semi-Inclusive Deep Inelastic Scattering (SIDIS). At variance with the hadronic collisions mentioned above, such process involves a large momentum transfer at the lepton vertex. The presence of a hard scale is a basic requirement in the derivation of a dedicated factorisation theorem [4,5] which ensures that collinear QCD factorisation holds in the leading-twist approximation for forward particle production in DIS. The relevant cross sections can then be factorised into perturbatively calculable short-distance cross sections and new distributions, fracture functions, which simulatenously encode informations both on the interacting parton and on the non-perturbative QCD dynamics of the spectator fragmentation into the observed forward hadron. Despite of being non-perturbative in nature, their scale dependence can be calculated within perturbative QCD [6]. Fracture functions obey in fact DGLAP [7] inhomogeneous evolution equations which result from the structure of collinear singularities in the target-fragmentation region [6,8]. The dedicated factorisation theorem [4,5] guarantees that fracture functions are universal distributions, at least in the context of SIDIS. On this theoretical basis, an impressive experimental program has been pursued at HERA in diffractive DIS which led to accurate determination of the so called diffractive parton distributions, i.e. proton-to-proton fracture functions in the very forward limit, allowing for the first time a quite accurate investigation of the parton content of the pomeron. The whole formalism has been later used in Ref. [9] to extract proton-to-Lambda fracture functions within a combined perturbative QCD fit to available SIDIS data.
In this paper we will focus on forward neutron production in DIS which provides, with respect to the aforementioned processes, complementary informations on soft QCD dynam- ics. An intensive physics program with forward neutron tagging has been pursued at HERA as well, where recent results [10,11] show that around 8% of the DIS events contain a forward neutron. These data are crucial in testing the limiting fragmentation hypothesis [12] and have been used as benchmark in a number of Regge-based models [13][14][15][16][17][18][19] which mainly concentrate on the modelisation of forward neutron production mechanisms. In the present paper we adopt instead a complementary approach and no modelisation of neutron production mechanisms is attempted. This strategy is in line with the factorisation theorem. The resulting set of proton-to-neutron fracture functions (nFFs) could then be used in hard-scattering factorisation test [20] in forward neutron tagged dijet photoproduction in ep collisions, as already measured at HERA [21,22], where factorisation is expected to hold only for the so-called direct component of the cross section. Even more intriguing appears to be the possibility of using nFFs for predicting the cross section for the associated production of a forward neutron and a Drell-Yan pair or dijet system in hadronic collisions. For this processes the factorisation theorem is not expected to hold and therefore nFFs determined by DIS data alone offer the opportunity to gauge factorisation breaking effects.
The paper is organised as follows. In Sec. II we describe the process under study, the kinematic variables relevant to the analysis and the observable used in the fit. In Sec. III we discuss the evolution of fracture functions and the general method with which we build initial conditions for the QCD evolution. In Sec. IV we describe the details of the QCD fit and in Sec. V we assess the impact of experimental and theoretical errors on the obtained neutron FF set. In Sec. VI we summarise our results.

II. DATA SET AND OBSERVABLE
In this analysis we consider semi-inclusive DIS events of the type where, beside the outgoing lepton, an additional neutron n is detected in the final state. In eq. (1) X stands for the unobserved part of the hadronic system and particles four-momenta are indicated in parenthesis. The kinematic variables Q 2 , x B and y are used to describe the inclusive DIS scattering process. They are defined as The kinematic variables used to describe the final state neutron are the neutron transverse momentum p T evaluated with respect to the beam axis and the longitudinal momentum fraction x L defined by where E n and E p are the neutron and proton energy in the laboratory frame, respectively.
In the following we use the scaled fractional momentum variable β defined by , respectively. In the one-photon exchange approximation, it reads:

III. THEORY SETUP
Hard-scattering factorisation for this class of processes states that the structure functions in eq. (5) are of the form The index i runs on the flavour of the interacting parton. The hard-scattering coefficients C ki (k = 2, L) are pertubatively calculable as a power expansion in the strong coupling α s and depend upon µ 2 F and µ 2 R , the factorisation and renormalisation scales, respectively. The C ki coefficient functions are the same as in fully inclusive DIS. The proton-to-neutron fracture functions M N i/P (β, µ 2 F , x L , p 2 T ) can be interpreted as the number density of interacting partons at a scale µ 2 F and fractional momentum β conditional to the observation of a forward neutron in the final state specified by a fractional momentum x L and transverse momentum squared T . They contain non-perturbative informations on the fragmentation of the spectator system which results from the hard interactions. The p T -unintegrated nFFs appearing in eq. (6) obey standard DGLAP [7] evolution equations [24]. In the case p T is integrated over up to values of order Q 2 , neutron fracture functions obey an inhomogenous DGLAP-type evolution equations [6]. The additional inhomogeneous term accounts for neutron production coming from the fragmentation of initial state parton radiation. In the present case, where the p T of the neutron is integrated up to some p T,max which lies in the non-perturbative region, neutron fracture functions are defined as and again obey familiar DGLAP evolution equations [25] valid at fixed values of x L . The central problem of this type of analyses is to find sensible initial conditions for the relevant distributions prior to evolution. One may resort to phenomenological models to describe forward neutron production. In general at low and intermediate x L the dominant mechanism is expected to be proton-remnant fragmentation into neutrons, while at high x L the exchange of virtual particles is expected to dominate. In the present analysis, we work at fixed x L and no attempts to model this non-perturbative dynamic at the proton vertex is made. Since hard-scattering factorisation in the form eq. (6) holds at fixed values of x L and p 2 T and this dependence is fully contained in nFFs, these conditional parton distributions are uniquely fixed by the kinematics of the outgoing neutron and they are, at least in principle, different for different values of x L and p 2 T . The approach we describe in this paper fully takes into account these important recipes in the construction of sensible input for nFFs distributions focusing on parton dynamics as explored by the virtual photon once an additional forward neutron is detected in the final state. This new approach have already been used in the extraction of diffractive parton densities from diffractive DIS data in Ref. [26]. This idea is realised in practice performing a series of QCD fits at fixed values of x L with a common initial condition controlled by a set of parameters {p i }. This procedure guide us to infere the approximate dependence of parameters {p i } on x L allowing the construction of a generalised initial condition in the (β, x L )-space to be used in a x L -combined QCD fit. It is important to note that if four-differential cross sections were available, the same method could be used to test whether, at fixed x L , the parton content explored by the virtual photon is the same (a part normalisation) in different neutron p T ranges.

IV. FITTING PROCEDURE
In this section we describe QCD fits at fixed values of x L . The distributions of neutron FFs in the quark sector at large β may show valence-like structures for some quark-flavour combinations. However the accessible values of β in the experimental data are quite low.
In view of this fact, and in order to reduce the number of free parameters, we assume that all light quark distributions are equal to each other, so that only the singlet and gluon distributions are required. We assume for the latter, at the arbitrary scale Q 2 0 and for any given value of x L , a momentum distributions of the type: structure functions at next-to-leading order which are then used to calculate the reduced cross section in eq. (5) and minimised against data by using the MINUIT program [28]. We adopt the generalised χ 2 definition proposed in Ref. [29] where systematics effects are incorporated in theory model predictions Here m i is the measurement of data point i, t i is the model prediction depending on a set of parameters p, σ i are the uncorrelated and statistical errors on data point i added in quadrature and ∆ ik is the correlated systematic error from source k on the data i. The variables s k denote Gaussian random variables with zero mean and unit variance. In the present Section we use the above definitions with s k = 0.
No cut on the invariant mass of the hadronic system X nor on the minimum Q 2 of data to be included in the fit is applied. As already discussed above, given the kinematic coverage of the data, we found that the singlet large-β coefficient C q is loosely constrained by data.
For the gluon distribution, which is only indirectly constrained by scaling violations, we found that B g strongly correlates with A g , when the former is left free to vary in the fit. For these reasons we set temptatively C q = 0.5, C g = 1 and B g = 0 so that the initial condition contains three free parameters in each x L -bin. We performed a combined scan on the value of initial scale Q 2 0 . Given the quite stiff functional form of the initial conditions, there is a mild dependence of the χ 2 on Q 2 0 which is then fixed at Q 2 0 = 1 GeV 2 . An essential condition for the x L -combination procedure to work is that good quality fits must be obtained in each x L -bin with the common initial conditions, eqs. (9,10). The χ 2 values of the fixed-x L fits, obtained with statistical and uncorrelated uncertainties added in quadrature, are presented in Tab. (I). From these values we may conclude that initial conditions provided by eqs. (9,10), supplemented by the constraints on C q , C g and B g , are general enough to describe the data in all x L -bins. With these results at hand we may now proceed and discuss the x L -combined fit. The generalised initial conditions have now the form where the β dependence is the same as in eqs. (9,10) and x L dependence is accounted for by the coefficients. The dependences of the A q , B q and A g free parameters on x L may be inferred inspecting Fig. (1). We adopt a redundant parametrisation of the coefficients of the where the term 1+c x d L is included to describe the relative maximum of the normalisation A q and A g at intermediate values of x L . For B q we assumed a second order polynomial in x L .
The parameters which were fixed in the fixed-x L fits are kept fixed to same values. All the   QCD settings are the same as in fixed-x L fits. With the help of eq. (13) and eqs. (14,15,16) we perform a series of x L -combined fits. At each iteration we study the eigenvalues of the covariance matrix of the fit parameters. Small eigenvalues, in fact, are associated to (combination of) parameters which are loosely determined by data. We found that, at large x L , both the singlet and gluon normalisations can be described by a common parametrisation The B q coefficient is found to be compatible with a constant so that only the parameter a 2 is left free to vary in the fit. Finally we found that b 1 is determined with rather large error and compatible with zero, so that we fix it to this value. The final form of the parametrisations of the coefficients is then for a total of seven free parameters. The best fit, with statistical and uncorrelated errors added in quadrature, returns a value of χ 2 =143. 6    approximation.

V. ERROR ESTIMATION AND PROPAGATION
In order to judge the agreement with other data sets and observables or to assess effects beyond the ones taken into account by the theoretical model, the obtained nFF parametrisation must be supplemented, fully exploiting the potential of the data, by a careful error analysis.
The general method with which experimental and theoretical uncertainties are propagated to a generic observable F is based on the construction of alternative nFFs parametrisation sets S k . By defining the difference r k = F (S k ) − F (S 0 ) and indicating with S 0 the best fit parametrisation, the uncertainties on a F are given by where θ is the Heaviside step function. In order to propagate statistical and uncorrelated uncertainties, following Refs. [30,31], we have diagonalised the covariance matrix of the best fit parameters and constructed a set of alternative parametrisations S k=1..14 according to the standard ∆χ 2 = 1 criterion. The error band constructed with the help of eqs. (21) and the S k=1..14 parametrisation set is shown in Fig. (1). The latter is narrower than individual errors on parameters obtained from fixed-x L fits. This error reduction is in fact due to the x L -combination, and can be understood considering for example B q , which is just a constant as a function of x L . In the combined fit, this parameter is constrained by 203 points rather than 29 of a single x L bin and so it is determined far more precisely. Since,however, fixed-x L fits, by construction, represent the best parametrisations of the data and eqs. (18,19,20) are interpreted as a mere interpolating tool, we require that the accuracy of the x L -combined fit does not exceeds the one of the fixed-x L fits. We found that a conservative ∆χ 2 = 9 criterion matches these requirements, as shown in Fig. (1). We also note that, incidentally, this number is close to the combination penalty reported in Tab. (II).
We now turn to the inclusion of systematics in the error analysis. In data from Ref. [23] nine systematics sources are identified plus the luminosity uncertainty [32]. For each of them we performed, according to the so-called offset method, alternative fits in which each s k is held fixed in turn either to -1 or +1 and produced the parametrisation set S k=15.. 34 . For some sources, for example the 5% luminosity uncertainty common to all data points, the shifts induce steady variation of the χ 2 and mostly correlates with the central values of the normalisations coefficients a 1 and a 3 . The impact of the propagation of systematic errors are presented in Fig. (3) where the best fit predictions are compared to ZEUS data [33].  [33]. The error bars associated with the data points show the sum in quadrature of the statistical and total systematic uncertainty. The light red error band corresponds to ∆χ 2 = 9 and it constructed with S k=1..14 . The grey error band is constructed as the quadratic sum of statistical (S k=1..14 ) and systematic (S k=15..34 ) contributions.
shown by the light grey band in Fig. (3). After taking into account all error sources we find that predictions based on nFFs describe ZEUS data both in shape and normalisation.
We conclude this section attempting an estimate of theoretical errors. Among them we mention the ones related to pQCD settings and the ones related to the choice of the parametrisation of initial conditions. The latter are by far dominant in the present analysis to the ∆χ 2 = 9 criterion, is obtained with the S k=1..14 and S k=35..38 parametrisation sets.
in the fit, i.e. C q = 0.5 and C g = 1. This implies that the error propagation produces, as shown in Fig. (4), an artificial shrinkage of the (light red) error band at large β and by no means represents a correct error estimate in this limit. In order to quantify the errors introduced by these assumptions, we performed four additional fits in which the C q and C g parameters are kept fixed in the minimisation but choosen in the following combinations: (C q , C g ) = (0, 1), (2, 1), (0.5, 0), (0.5, 3). The latter values are choosen such that fits return a difference of around ∆χ 2 = 9 with respect to best fit. The latter four alternative parametrisations, S k=35..38 , are used to produce the light grey errorband shown in Fig. (4) and shows the degree of underdetermination of the distributions in this region.
These additional parametrisations can be especially useful to propagate uncertainties to observables which require nFFs large-β extrapolation, for example jet cross section. Less problematic appear the Q 2 extrapolation, since the latter is fully predicted by the theory.
In Fig. (5) we present the initial condition at x L = 0.635 for three different values of Q 2 . It is interesting to note that the uncertanties on the nFFs increase as Q 2 decrease. This effect can be partly ascribed to the fact that no data point with Q 2 below 7.3 GeV 2 is included in the fit. But, more importantly, it has to be ascribed to QCD evolution: small displacements of the parametrisations at high Q 2 , where they are actually constrained by the data, turn into large fluctuations of the initial conditions at Q 2 0 , due to the logarithmic nature of QCD evolution equations.

VI. CONCLUSIONS
In this paper we have presented a perturbative QCD analysis and extraction of neutron fracture functions from forward neutron production in DIS in HERA kinematics. Data can be decribed by the leading-twist approximation implied by the hard-scattering factorisation formula and perturbative QCD evolution down to the lowest values of Q 2 accessed by the experiment. The results of the fit and within the precision of the present data, indicate that the proton-vertex factorisation hypothesis is supported to good accuracy, a fact which is likely to be related to the relative low β regime accessed by the measuraments. The nFFs low-Q 2 extrapolation, although with large uncertainties, can be used to address the impact of absorptive effects going from the DIS to the photoproduction regime. The predictions based on the obtained nFFs have been succesfully compared to ZEUS data and an error estimation on nFFs has been provided. The obtained nFFs parametrisation obtained in NLO QCD is a quantitative tools which can be used in factorisation tests in processes with a tagged forward neutron. In this context we mention dijet photoproduction in ep collisions and dijet or Drell-Yan pair production in hadronic collisions. For the latter process, next-toleading order corrections have been estimated in Ref. [34]. The nFFs parametrisation and error set are available upon request to the author and are provided as fortran steering file for the QCDNUM17 package.