Impact of the LHCb $p\!+\!\!^{4}\!H\!e$ fixed-target $D^0/\bar{D^0}$ data on the intrinisic $c \bar c$ component in the nucleon

We discuss the impact of the new LHCb fixed-target $p\!+\!^{4}\!He$ data for $D^0/{\bar D}^0$ production on the intrinsic $c \bar c$ component in the nucleon wave function. Within the scenario presented here neither the traditional gluon-gluon fusion or quark-antiquark annihilation mechanisms calculated in the $k_T$-factorization approach nor their counterparts from the collinear next-to-leading order collinear framework are sufficient to describe the transverse momentum and rapidity distributions of $D^0/{\bar D}^0$ mesons. First the $c{\bar c}$-pair production within the standard frameworks is considered. Here a crucial role of the $c \to D$ hadronization effects at low energies and low transverse momenta is found and discussed, which was not analyzed in previous studies. A contribution related to intrinsic $c \bar c$ component in the nucleon wave function is included in addition. Two models of the symmetric ($c(x) = \bar c(x)$) intrinsic charm (IC) component are considered. The intrinsic charm $g^* c \to g c$ (or $g^* {\bar c} \to g {\bar c}$) contribution needs to be regularized in order to obtain a suppression of the minijet $p_{T}$ spectrum present in the phenomenological minijet model, commonly used in Monte Carlo generators. We show that in our model the regularization parameter can be obtained from the fit to the LHCb fixed-target data under consideration here. We discuss uncertainties of our calculations (scale, charm quark mass, fragmentation function) as well as set limits on the IC probability. According to our model the intrinsic charm probability $P_{IC} = 1.65\%$ allows to significantly improve description of the LHCb data but the number is rather uncertain.


I. INTRODUCTION
The proton is known to be a complicated object built of QCD degrees of freedom: quarks, antiquarks and gluons or even their conglomerates. However, its wave function remains to large extent unknown. The flavour structure of the nucleon is especially interesting.
An interestingd −ū asymmetry confirmed by the Gottfried sum rule violation [1] and fixed target Drell-Yan experiment [2] is one example. This is usually explained in terms of the pion (meson) cloud in the nucleon (see e.g. Ref. [3]). The pion cloud probability is of the order of 30%. It contributes to nonperturbative light quark/antiquark distributions in the nucleon. Another example is the strangeness in the nucleon. Here the kaon cloud leads to a nonperturbative contribution of s/s quarks/antiquarks. A few percent effect of the nonperturbative ss component was predicted [3].
The situation for the charm-anticharm content of the nucleon is less understood. Long ago Brodsky et al. [4] proposed a simple model of the large-x intrinsic charm due to uudcc Fock component of the wave function. As a consequence of heavy mass of the c-quark orc-antiquark their longitudinal momentum fractions are of the order of x ∼ 0.3 − 0.5. In its simplest form the model leads to c(x) =c(x) (symmetric IC). Another possible mechanism is nonperturbative gluon splitting into cc. Such a mechanism was discussed e.g. in Refs. [5,6]. Also here c(x) =c(x), but typical longitudinal momentum fractions of c-quark orc-antiquark are much smaller than in the Brodsky-Hoyer-Peterson-Sakai (BHPS) model. We shall call the second model as the SEA-like. It is very difficult to calculate in a realistic way the integral c(x)dx = c(x)dx (probability of such a component in the wave function) as it goes beyond the scope of perturbative methods.
In Ref. [7] the c(x) andc(x) instrinsic charm/anticharm distributions are calculated postulating heavy charmed meson -heavy charmed baryon component in the nucleon wave function. In this approach typically c(x) is somewhat different thanc(x) (asymmetric IC) but the most probable range of x is similar as for the BHPS model. Therefore our analysis of the BHPS model is approximately valid also for the meson cloud approach.
The intrinsic charm component was discussed recently in Refs. [8][9][10] in the context of the associated production of prompt photons and charm-quark jets in pp-collisions at LHC energies. In Refs. [11][12][13] an associated production of Z 0 and heavy flavoured (HF) jets was studied in the context of identifying intrinsic charm in the nucleon. Recently, in Ref. [13] the authors found that the forthcoming ATLAS and CMS measurements of Z 0 + HF production at √ s = 13 TeV can be very important in the context of searching for the IC contribution in the proton. They suggested measuring a new observable, defined as the double ratio of cross sections σ f orward (Z+c) (where central and forward means |y Z | < 1.5 and 1.5 < |y Z | < 2.5, respectively). This was found to be extremely sensitive to the IC signal, which would be potentially visible at large transverse momenta of Z 0 -boson and/or heavy flavour jet. At low energy discussed in the present paper the production of Z 0 is not possible. On the other hand registration of the extra photon is not possible due to very low luminosity in the fixed target experiment.
In the present study we are interested rather in production of c-quark and/orcantiquark and associated D mesons in the forward direction. The interaction of gluon (small-x) from the projectile with c orc (large-x) in the target ( 4 He) is the underlying partonic mechanism [6]. The 4 He is rather light nucleus and the associated nuclear effects are of the order of a few percent only and can be safely neglected. We shall discuss this issue shortly in the main text. Focusing on rather small charm transverse momenta, the cg * → cg (orcg * →cg) subprocess needs to be regularized in order to obtain a phenomenologically motivated suppression of the minijet p T spectrum. Here we follow the minijet model originally proposed in Ref. [14] and further adopted in modern Monte Carlo generators (see e.g. PYTHIA [15]). The parameter governing the regularization must be adjusted to experimental data. It will be discussed how to do it in our case.
In parallel we have made a similar analysis of the IC signal for high-energy prompt neutrino production in the Earth's atmosphere coming from the decay of charmed meson [16]. The two analyses (here and for atmospheric neutrino) are complementary. Combining conclusions from both of them may be therefore very useful.

A. The leading perturbative mechanism
We remind here very briefly the theoretical formalism for the calculation of the ccpair production in the k T -factorization approach [17], as adopted and discussed by one of us in the context of the LHCb fixed-target charm data in Ref. [18]. In this framework the transverse momenta k t 's (virtualities) of both partons entering the hard process are taken into account, both in the matrix elements and in the parton distribution functions.
Emission of the initial state partons is encoded in the transverse-momentum-dependent (unintegrated) PDFs (uPDFs). In the case of charm flavour production the parton-level cross section is usually calculated via the 2 → 2 leading-order g * g * → cc fusion mechanism with off-shell initial state gluons that is the dominant process at high energies (see Fig. 1). Even at lower energies as long as small transverse momenta and not extremely backward/forward rapidities are considered the q * q * → cc mechanism remains subleading. Then the hadron-level differential cross section for the cc-pair production, formally at leading-order, reads: where F g (x 1 , k 2 1,t , µ 2 F ) and F g (x 2 , k 2 2,t , µ 2 F ) are the gluon uPDFs for both colliding hadrons and M off−shell g * g * →cc is the off-shell matrix element for the hard subprocess. The gluon uPDF depends on gluon longitudinal momentum fraction x, transverse momentum squared k 2 t of the gluons entering the hard process, and in general also on a (factorization) scale of the hard process µ 2 F . The extra integration is over transverse momenta of the initial partons. Here, one keeps exact kinematics from the very beginning and additional hard dynamics coming from transverse momenta of incident partons. Explicit treatment of the transverse momenta makes the approach very efficient in studies of correlation observables. The two-dimensional Dirac delta function assures momentum conservation. The gluon uPDFs must be evaluated at longitudinal momentum frac- i,t + m 2 c is the quark/antiquark transverse mass. As we have carefully discussed in Ref. [19], there is a direct relation between a resummation present in uPDFs in the transverse momentum dependent factorization and a parton shower in the collinear framework. Such observations have been made also earlier, for example in Refs. [20,21]. The commonly accepted statement is that actually in the k T -factorization approach already at leading-order some part of radiative higherorder corrections can be effectively included via uPDFs, without any additional showering procedure. However, it is true only for those uPDF models in which extra emissions of soft and even hard partons are encoded, including k 2 t > µ 2 F configurations. In most uPDF cases the off-shell gluon can be produced either from gluon or quark, therefore, in the k T -factorization all channels driven by gg, qq and even by qg initial states are open already at leading-order (in contrast to the collinear factorization). Then, when calculating the charm production cross section via the g * g * → cc mechanism one could expect to effectively include contributions related to an additional one or two (or even more) extra partonic emissions which in some sense plays a role of the initial state parton shower.
The kinematical configuration of the fixed-target LHCb experiment corresponds to the region where in principle the Catani-Ciafaloni-Fiorani-Marchesini (CCFM) [22] evolution equation is legitimate for any pQCD theoretical calculations and could, in principle, be used to describe the dynamics behind the mechanisms of e.g. open charm meson production. Within the CCFM approach the parton transverse momentum is allowed to be larger than the scale µ 2 . This useful feature translates into the easy and effective taking into account of higher-order radiative corrections, that correspond to the initial-state real gluon emissions which are resummed into the uPDFs. In the numerical calculations below we follow the conclusions from Ref. [18] and apply the gluon uPDFs obtained from the CCFM evolution equation. There a different models of the CCFM unintegrated gluon densities have been tested and found to lead to a reasonable description of the LHCb fixed-target charm data. The best theory to data relation has been obtained in the case of the most up-to-date JH-2013-set1 and JH-2013-set2 unintegrated gluon densities [23] that are determined from high-precision DIS measurements.
As a default set in the numerical calculations we take the renormalization scale µ 2 = it n (averaged transverse mass of the given final state) and the charm quark mass m c = 1.5 GeV. The strong-coupling constant α s (µ 2 R ) at next-to-next-to-leading-order is taken from the CT14nnloIC PDF routines 1 [24]. The CCFM uPDFs here are taken at rather atypical value of the factorization scale µ 2 F = M 2 cc + P 2 T , where M cc and P T are the cc-invariant mass (or energy of the scattering subprocess) and the transverse momentum of cc-pair (or the incoming off-shell gluon pair). This unusual definition has to be applied as a consequence of the CCFM evolution algorithm [23]. For completeness of the present study, in some places we also apply the Martin-Ryskin-Watt (MRW) model for unintegrated gluon densities [25]. The intrinsic charm contribution to charm production cross section (see Fig. 2) is obtained within the hybrid theoretical model discussed by us in detail in Ref. [6]. The LHCb fixed-target configuration allows to explore the charm cross section in the backward rapidity direction where an asymmetric kinematical regime can be explored. Thus in the basic gc → gc reaction the gluon PDF and the intrinsic charm PDF are simultaneously probed at different longitudinal momentum fractions -rather intermediate for the gluon and large for the charm quark.

B. The intrinsic charm contribution
Within the asymmetric kinematic situation x 1 ≪ x 2 the cross section for the processes 1 In the numerical calculations in Ref. [18] the strong-coupling at leading-order was taken. Here, in order to consistently compare the standard mechanism and the intrinsic charm component we keep in both cases the α s (µ 2 R ) values encoded in the CT14nnloIC PDF.
under consideration can be calculated in the so-called hybrid factorization model motivated by the work in Ref. [27]. In this framework the small-or intermediate-x gluon is taken to be off mass shell and the differential cross section e.g. for pp → gcX via g * c → gc mechanism reads: is the unintegrated gluon distribution in one proton and c(x 2 , µ 2 ) a collinear PDF in the second one. The dσ g * c→gc is the hard partonic cross section obtained from a gauge invariant tree-level off-shell amplitude. A derivation of the hybrid factorization from the dilute limit of the Color Glass Condensate approach can be found in Ref. [28] (see also Ref. [29]). The relevant cross sections are calculated with the help of the KATIE Monte Carlo generator [30]. There the initial state quarks (including heavy quarks) can be treated as a massless partons only.
Working with minijets (jets with transverse momentum of the order of a few GeV) requires a phenomenologically motivated regularization of the cross sections. Here we follow the minijet model [14] adopted e.g. in PYTHIA Monte Carlo generator, where a special suppression factor is introduced at the cross section level [15]: for each of the outgoing massless partons with transverse momentum p t , where p T0 is a free parameter of the form factor that also enters as an argument of the strong coupling constant α S (p 2 T0 + µ 2 R ). This suppression factor was originally proposed to remove singularity of minijet cross sections in the collinear approach at leading-order. In the hybrid model (or in the k Tfactorization) the leading-order cross sections are finite as long as k T > 0, where k T is the transverse momentum of the incident off-shell parton. Within this approach, a treatment of the small-k T region in the construction of a given unintegrated parton density is crucial.
Different models of uPDFs may lead to different behaviour of the cross section at small minijet transverse momenta but in any case the cross sections should be finite. However, as it was shown in Ref. [31], the internal k T cannot give a minijet suppression consistent with the minijet model and the related regularization seems to be necessary even in this framework.
In the numerical calculations below, the intrinsic charm PDFs are taken at the initial scale m c = 1.3 GeV, so the perturbative charm contribution is intentionally not taken into account. In the numerical calculations below we apply different grids of the intrinsic charm distributions from the CT14nnloIC PDF [24] that correspond to the BHPS model [4] as well as to the so-called SEA-like models. In Fig. 3 we compare the collinear gluon (left panels) ant the charm quark (right panels) PDFs obtained with (CT14nnloIC) and without (CT14nnlo) the intrinsic charm concept at different factorization scales (upper and lower panels). Here, the intrinsic charm grid that correspond to the BHPS model with 1% probability for intrinsic charm was used. In this case the intrinsic charm contribution leads to a significant enhancement of the charm PDF in the region of longitudinal momentum fraction x > 0.1.

C. The charm quark to meson transition
The transition of charm quarks to open charm mesons is done in the framework of the independent parton fragmentation picture (see e.g. Refs. [32,33]). In Ref. [18] we followed the standard prescription, where the inclusive distributions of open charm meson were obtained through a convolution of inclusive distributions of charm quarks/antiquarks and c → D fragmentation functions: where p t,c = p t,D z and z ∈ (0, 1). There the typical approximation was done that y c is unchanged in the fragmentation process, i.e. y D = y c . This commonly accepted and frequently used method was originally proposed for light partons. It can be safely used only when both, mass of the parton and mass of the hadron can be neglected [33].
In principle, this approximation may not be valid for the case of heavy and even light parton fragmentation to heavy object, especially, at lower energies or/and considering regions of small meson transverse momenta. It is obvious that working with massive particles this model may break down when approaching p T ∼ m D region. In this regime one could expect a violation of "energy conservation" and events with hadrons that have larger energies than the energy of the parent parton can frequently appear. In some corners of the phase space the E D < E c relation may be broken very strongly. As long as one is considering c.m.s. midrapidities and/or large c.m.s. collision energies this mass effect shall be rather negligible, especially, when a low transverse momentum cut is applied.
However, the situation may dramatically change when going to lower energies and discussing backward or forward production. This might be exactly the case of the LHCb fixed-target charm data.
Therefore, here we also follow a different idea and assume that the D-meson is emitted in the direction of parent c-quark/antiquark, i.e. η D = η c (the same pseudorapidities or polar angles). Within this approach still different options for z-scaling come into game, including e.g. the three-momentum (p c = p D z ) and the light-cone momentum (p [32].
In numerical calculations we take the Peterson fragmentation function [34] with ε = 0.05, often used in the context of hadronization of heavy flavours. Then, the hadronic cross section is normalized by the relevant charm fragmentation fractions for a given type of D meson [35]. In the numerical calculations below for c → D 0 meson transition we take the fragmentation probability P c→D = 61%.

III. NUMERICAL RESULTS
Let us start presentation of our numerical results with the differential cross sections for charm quark production in pp-scattering at √ s = 86.6 GeV. In Fig. 4 we show the transverse momentum (left panels) and rapidity (right panels) distributions of charm quark for the standard gg → cc partonic mechanism (top panels) as well as for the gc → gc mechanism with intrinsic charm in the initial state (bottom panels). Considering the standard charm production mechanism the results are obtained within the k T -factorization approach for different models of gluon uPDFs, namely the JH2013set1 (long-dashed histograms), the JH2013set2 (solid histograms), and the MRW-CT14lo (dashed histograms).
These results are compared in addition to the next-to-leading order collinear predictions obtained according to the FONLL [36] framework (dotted histograms). Numerical results for the mechanism driven by the intrinsic charm are obtained within the hybrid approach with off-shell gluon and with collinear intrinsic charm distribution taken from the CT14nnloIC PDFs. Here we use the PDF set that correspond to the BHPS model with 1% probability for finding intrinsic charm content in a proton. Here, the hybrid model results are compared with the leading-order collinear predictions. We observe that the different models of the gluon uPDF lead to quite different charm quark distributions. The statement is valid for both considered partonic mechanisms. Clearly, the unintegrated gluon densities are a source of large uncertainties of the predictions (see also discussion in Ref. [6]). Therefore their application in the probed kinematics needs an experimental verification. In the case of the standard charm production mechanism it has been already done very recently in Ref. [18].  In the present paper, we wish to concentrate mostly on the mechanism associated with the intrinsic charm concept in proton which was not discussed so far in the context of the LHCb fixed-target charm data. Before we go to the main results we wish to visualize qualitatively the kinematics behind the considered production mechanism. In Fig. 5 we present double differential parton-level cross section for charm quark as a function of longitudinal momentum fractions log 10 (x 1 ) and log 10 (x 2 ) carried by the incident partons for cg * → gc (left panel) and g * c → gc (right panel) mechanisms. Here we impose on  Here results for three different values of p T0 are shown. Details are specified in the figure.
As was already introduced in the previous section (and also discussed in Ref. [6]), our model for calculation of the intrinsic charm contribution to the charm production cross section depends on a free parameter p T0 used for a regularization of the cross section. From Fig. 6 we see that the predictions for charm quark transverse momentum Ref. [26] (see Fig. 4 therein). Fig. 8 shows the ratio of nuclear ( 4 He) to nucleon gluon distributions for a typical for our reaction, rather low, factorization scale µ 2 = 5 GeV 2 . The effect is rather small and different models give rather contradictory results. Therefore in the following we shall neglect the nuclear effects completely. Below we show in addition the theory/data ratio. We start with the results for the standard cc-pair production mechanism calculated both, in the framework of the k T -factorization approach with offshell initial state partons and with the CCFM uPDFs, and in the NLO collinear approach within the FONLL framework. In Fig. 7 we present differential distributions of D 0 meson as a function of c.m.s. rapidity (left panel) and transverse momentum (right panel) for p + 4 He collisions together with the LHCb experimental data points [26]. A detailed phenomenological study of the data set within the k T -factorization approach has been already performed by one of us in Ref. [18]. There different models of unintegrated gluon densities were used and tested against the LHCb data. It was found that, in principle, the calculations based on the k T -factorization with the CCFM unintegrated gluon densities are able to describe the data in a satisfactory manner. Here we repeat a part of the calculations presented there and use two different sets of the gluon CCFM uPDFs: the JH-2013-set1 and JH-2013-set2 fits [23], however, the results obtained here may slightly differ from those reported in Ref. [18] since the choice of the strong coupling is not the same. Here and in the following, we keep the α S (µ 2 ) at NNLO in order to be consistent with the α S (µ 2 ) grids present Now we wish to discuss the aspect of above calculations that was not analysed in Ref. [18]. Namely, we wish to pay more attention to the details of the fragmentation procedure. In the following we show for the first time the result of the η D = η c prescription with scaling variable defined as follows: p H = zp q (momentum scaling) or (E + p) H = z(E + p) q (light-cone scaling). In Fig. 9 we compare results obtained within these two prescriptions to those calculated within the standard procedure with unchanged rapidity. Clearly, a significant difference between the standard and the lowenergy model appears at small meson transverse momenta (p T < 2 GeV) which in consequence leads to a quite different rapidity distribution. It seems that in the backward rapidity region details of the fragmentation procedure are of a special importance and could lead to different conclusions about quality of the theoretical predictions and understanding of the experimental data.
As it is shown in Fig. 11, when applying the low-energy model for fragmentation with the light-cone momentum scaling the overall picture presented in Fig. 7 is changed. . Different fragmentation procedures lead to different scenarios. The new scenario presented here might be more accurate than the standard one in the kinematics probed by the LHCb in its fixed-target mode. The fact that within this model the LHCb data favour more results for the JH-2013-set1 gluon uPDF than for the JH-2013-set2 gluon uPDF is a bit unexpacted, since the latter is determined from the fit to both F (charm) 2 and F 2 DIS data.
Considering the FONLL results plotted in Fig. 11 (bottom panels) we see that the cen- tral predictions visibly underestimate both data sets when m c = 1.5 GeV is taken. The situation is slightly better for a smaller charm quark mass, i.e. m c = 1.3 GeV, however, this choice should be rather considered as an upper limit and not a central prediction.
Keeping m c = 1.5 GeV, we observe also that there is a sizeable difference between the FONLL prediction and the k T -factorization result in the very backward rapidity region.
We find the observation to be a direct consequence of the difference in the gluon densities used in both frameworks rather than the effect related to the subleading qq → cc production channels, part of which might be missing in the k T -factorization approach 2 .
As we checked numerically at the considered √ s = 86.6 GeV collision energy, for both the LO and the NLO collinear calculations the gg-fusion mechanism dominates over the qqannihilation channels even in the far forward/backward rapidity direction (see Fig. 10).
Starting from this point let us now follow a scenario in which the prediction obtained within the k T -factorization framework with the JH-2013-set2 gluon uPDF, including fragmentation with the light-cone momentum scaling, is the point of reference suggested by the theory behind the model ingredients. Potentially, this picture leaves a room for 2 In principle, according to the qq → cc production mechanism, contributions coming from the light valence quark interactions might be limited in the k T -factorization framework (in contrast to the sea-quark interactions). Before we start a presentation of the predicted intrinsic charm contribution to the differential distributions reported by the LHCb experiment we shall check what is a contribution to the charm production cross section coming from the g * c → gc mechanism under the assumption that there is no non-perturbative intrinsic cc content in the nucleon. So we test here the impact on the charm cross section of the charm quarks in the initial state produced only perturbatively, during the QCD evolution. Here we also take the charm quark PDF at the fixed factorization scale µ = 1.3 GeV so we keep the PDF not much evolved from the initial condition. The corresponding results of our calculation are shown in Fig. 12. Clearly, the g * c → gc mechanism with no intrinsic charm leads to a smaller cross section than the one predicted by the standard (gluon-gluon fusion or quark-antiquark annihilation) cc-pair production mechanism and cannot improve description of the LHCb open charm fixed-target data.
In the next step we examined also contributions coming from interactions of quarks only. As was shown e.g. in Ref. [13], the mechanisms driven by qq and cq interactions are competitive with respect to the gluon-gluon fusion in the case of associated production of Z-boson and heavy flavoured jets at high energies. These contributions were found there to be important in the context of intrinsic charm studies. Here we present a similar analysis and include into numerical calculations the qq * → cc (dash-dotted histograms) and the cq * → qc (dashed and dotted histograms) mechanisms. The latter is calculated in two different ways: when only non-perturbative intrinsic charm quark PDF at the initial scale is taken into account (dotted histograms) as well as when only the perturbative  charm quarks generated in the DGLAP evolution are used (dashed histograms). Here the PB-NLO-set1 [37] unintegrated quark distributions are used for the small-x region.
As can be concluded from Fig. 13, at the energy and kinematics considered here, all three contributions are found to be subleading with respect to mechanisms with gluons in the initial state and can be safely neglected in the present study. Even in the very backward rapidity bin the quark initiated contributions are order of magnitude smaller than the leading mechanisms with gluons in the initial state.
Let us now show how the situation changes when the intrinsic charm content of the nucleon is taken into account. We start with predictions obtained within the SEA-like model for the intrinsic charm. In Fig. 14   transverse momentum rather than at larger p T -valus. The prediction of the SEA-like P = 1.6% overshoots the LHCb data significantly at small meson transverse momentum.
It also overestimates the rapidity spectrum in the whole considered range. The SEA-like P = 0.6% result lies closer to the data points, giving a reasonable description but still seems to overestimate the low-p T region and do not improve description of the data at the large-p T part of the distribution. As a next step we took into consideration the BHPS model of the intrinsic charm. In to the BHPS model with the P = 0.6%, P = 1.0%, and P = 2.1%, respectively. In contrast to the SEA-like case, the predictions of the BHPS model contribute significantly to the considered cross section especially at larger meson transverse momentum and in the region of the backward rapidity. Here the intrinsic charm contribution does not really affect the standard results at midrapidity and small meson transverse momentum which is more supported by the LHCb experimental data. A presence of the intrinsic charm component improves description of the data within the considered scenario. We see that the P = 0.6% scenario is favoured by the LHCb fixed-target data. However, as can be seen from Fig. 16, this conclusion depends on the model parameter p T0 used in the calculation. For the larger value of the regularization parameter, i.e. p T0 = 2.0 GeV, the P IC = 1.0% scenario seems to be more appropriate. This sensitivity almost disappears for large meson transverse momentum. The predicted intrinsic charm contribution for the last bin in p T does not really depend on the choice of the p T0 what is explicitly shown in Fig. 17. Therefore the last point of the measured p T spectrum can be used to constrain the P IC probability. On the other hand, having the P IC fixed by the data for transverse momentum distribution, the corresponding spectrum in rapidity (especially the very last backward point) could be used to constrain the regularization parameter p T0 , which is the only free parameter in our model. The rapidity distribution measured by the LHCb experiment strongly supports larger values of the regularization parameter, namely the p T0 = 2.0 GeV.
As one can see from above, each of the three different values of the intrinsic charm probability visibly modify the p T -slope of the differential distribution as well as the backward rapidity region. A presence of the IC components improves quality of the experimental data description. Within the presented theoretical model it is possible to roughly set the upper limit for the intrinsic charm probability to be ≈ 1% which is consistent with the central prediction of the CT14nnloIC PDF global fit.  We wish to discuss now uncertainties related to our calculations. Below we discuss the scale, mass (Fig. 18) and the fragmentation function (Fig. 19) uncertainties as well as the χ 2 analysis (Fig. 20) used in order to set limits on the IC probability.
In Fig. 18 we show uncertainties related to quark mass and the choice of the scale.
The uncertainties are not small, however the shapes of distributions, particularly that for dσ/dy, are quite different than the experimental ones and the IC contribution visibly improves description of the rapidity distribution.
In Fig. 19 we show the uncertainties due to the choice of the fragmentation functions.
Such uncertainties are much smaller than those due to scale or charm quark mass and can be safely neglected.
In Fig. 20 we show χ 2 as a function of the intrinsic charm probability P IC . Here we follow the method presented in Refs. [8,9] and define χ 2 as where y i is the measurement, f (P IC ) i is the theoretical result for a given P IC and ∆σ i is is the sum in quadrature of the uncertainties coming from the measurement and the  uncertainties coming from the theoretical calculations. We observe a minimum at P IC ∼ 1.65 % but the uncertainty is really large. A large range of the IC probability is still allowed. In parallel we made a study of extracting the IC component from the IceCube neutrino data. There only upper limit P IC 2.0% could be extracted [16]. This seems consistent with the present result. Similar result was obtained in Ref. [8].

IV. CONCLUSIONS
In the present paper we have discussed in detail how the relatively recent data on D 0 /D 0 production in p+ 4 He fixed-target LHCb experiment put constraints on intrinsic charm component in the nucleon. We have considered two models of cc symmetric intrinsic charm: the well known BHPS model as well as a SEA-like model. In our approach the knocked out c-quark orc-antiquark fragment producing D 0 orD 0 mesons. This part of the reaction is described by a phenomenological Peterson fragmentation function correctly using kinematics.
As a "background" contribution we have used the traditional gluon-gluon and quarkantiquark annihilation mechanism. The respective contributions were calculated within the k T -factorization approach. The gluon-gluon fusion was calculated using so-called CCFM unintegrated gluon distribution adjusted to experimental F 2 structure function.
The quark-antiquark annihilation contributions turned out negligible compared to the sum of gluon-gluon and IC contributions.
Within the scenario presented here the traditional components seem insufficient to describe the LHCb data, especially when discussing the very backward rapidity region and large meson transverse momenta. Here we follow a different treatment of the fragmentation process with respect to the predictions presented by one of us in Ref. [18]. A corresponding consequences related to the modelling o hadronization effects at lower energies are carefully discussed.
We have studied uncertainties of our calculation related to scale and charm quark mass as well as c → D fragmentation functions. The first two are large. However, the shape of rapidity distribution of the traditional calculation is rather different than the experimental one. The uncertainty due to fragmentation functions seems rather small. The BHPS IC improves the description of the data. We have also tried to set limits on the probability of the IC component. We have observed a minimum of χ 2 at P IC ≈ 1.65 %. There is, however, a big uncertainty and the probability 0.2-3.15 % within 1σ confidence level.
For comparison the upper limit obtained from the analysis of the neutrino IceCube data [16] is P IC 2.0% , which is consistent with the present analysis. Similarly, in Refs. [8,9] the authors found P IC 1.93% using tha ATLAS data on the associated production of prompt photons and charm-quark jets.