$J/\psi$ meson production in association with an open charm hadron at the LHC: a reappraisal

We critically (re)examine the associated production process of a $J/\psi$ meson plus an open charm hadron at the LHC in the proton-proton ($pp$) and proton-lead ($p{\rm Pb}$) collisions. Such a process is very intriguing in the sense of tailoring to explore the double parton structure of nucleons and to determine the geometry of partons in nuclei. In order to interpret the existing $pp$ data with the LHCb detector at the centre-of-mass energy $\sqrt{s}=7$ TeV, we introduce two overlooked mechanisms for the double parton scattering (DPS) and single parton scattering (SPS) processes. Besides the conventional DPS mode, where the two mesons are produced almost independently in the two separate scattering subprocesses, we propose a novel DPS mechanism that the two constituent (heavy) quarks stemming from two hard scatterings can form into a composite particle, like the $J/\psi$ meson, during the hadronisation phase. However, it turns out the corresponding contribution is small in $J/\psi+c\bar{c}$ hadroproduction. On the contrary, we point out that the resummation of the initial state logarithms due to gluon splitting into a charm quark pair is crucial to understand the LHCb measurement, which was overlooked in the literature. We perform a proper matching between the perturbative calculations in different initial-quark flavour number schemes, generically referring to the variable flavour number scheme. The new variable flavour number scheme calculation for the process strongly enhance the SPS cross sections, almost closing the discrepancies between theory and experiment. Finally, we present our predictions for the forthcoming LHCb measurement in $p{\rm Pb}$ collisions at $\sqrt{s_{NN}}=8.16$ TeV. Some interesting observables are exploited to set up the control regions of the DPS signal and to probe the impact-parameter dependent parton densities in lead.


I. INTRODUCTION
The associated production of a J/ψ meson plus an open charm hadron at the LHC is an intriguing process. It provides a mean to carry out studies on several QCD phenomena. First of all, the very large yields of both charm (anti-)quark and charmonium give rise naturally to the emergence of two simultaneous parton scattering subprocesses, aka double parton scattering (DPS), in a single proton-proton collision due to the composite nature of the proton. In the conventional single parton scattering (SPS), the process has also been proposed to study the heavy quarkonium production mechanism [1][2][3][4][5][6], the intrinsic charm content of the proton [7,8] and the factorisation-breaking effect due to colour transfer between the J/ψ and the open charm in the mass threshold regime [9,10]. The measurement of such a process in heavy-ion collisions is further motivated to understand the spatial dependence of nuclear modification of parton flux in nuclei [11], given the observation of the large nuclear modification in the inclusive productions of J/ψ and open charm (see, e.g., Ref. [12] and references therein).
Despite the numerous theoretical motivations, the only measurement at the LHC so far was carried out with the LHCb detector based on a relatively low integrated luminosity 355 pb −1 of Run-I proton-proton (pp) data at √ s = 7 TeV in 2012 [13]. Yet, the measurement is precise enough to be conclusive in the yields. It is claimed that the SPS stems from the gluon fusion gg → J/ψ + cc (thereafter referred to as "ggF") is quite negligible compared to the observed production rate. For instance, the measured cross section of J/ψ + D 0 and its charge conjugate mode J/ψ +D 0 in the LHCb acceptance [2.0 < y J/ψ , y C < 4.0, P T,J/ψ < 12 GeV and 3 < p T,C /GeV < 12, where C is either a charm hadron (D 0 , D + , D + s , Λ + c ) or its charge conjugate] amounts to 161.0 ± 3.7 ± 12.2 nb, while the corresponding ggF cross section is significantly smaller ranging from 3.7 to 16 nb. However, the DPS dominance picture seems to be in trouble in interpreting the shapes of the differential distributions, in particular for the transverse momentum P T of J/ψ, by assuming the standard zero correlation hypothesis between the J/ψ and D 0 mesons. Under the simple assumption, the DPS cross section amounts to where we have denoted DPS as DPS 1 in order to differentiate another DPS mechanism that will be introduced in Sec. II. Presumably, by construction, the shape of such a DPS 1 distribution should be identical to that of the single inclusive (prompt) J/ψ production, while the observed spectrum is much harder than the corresponding DPS 1 prediction (see Fig.11a in Ref. [13]). 1 In fact, the measured P T,J/ψ spectrum matches the ggF SPS result simulated with HELAC-Onia 2.0 [14,15] and Pythia 8.186 [16] as reported in Fig. 1a. We have used the normalisation of DPS 1 with σ eff,pp = 15 mb for the sum of both contributions (black solid line), which correctly reproduces the total yield, and the band stands for the missing higher-order QCD radiative corrections using the standard scale variation. Similar observation can be found in the invariant mass distribution of the two meson system in Fig. 1b. While there is no obvious reason that the zero correlation assumption in DPS 1 is strongly violated by the initial gluon-gluon correlation 2 since this assumption works pretty well in other processes except the normalisation encoded in σ eff,pp , it is clear that a coherent physical picture is still missing in interpreting these LHCb data, and therefore prevents us from understanding the process and from using the process as a tool to study other interesting phenomenology.
The main purpose of the present paper is to address the above mentioned issue in pp collisions and to present our predictions for the forthcoming LHCb measurement of the same process in proton-lead (pPb) collisions at √ s N N = 8.16 TeV. In addition, we would like to deduce the potential of constraining the impact-parameter dependent parton densities in nuclei by measuring the process at the LHC. The remaining context is structured as follows. We introduce a new DPS mechanism in Sec. II, and discuss the theoretical aspects of the initial (anti-)charm contributions and the variable flavour number scheme to improve our SPS predictions in Sec. III. In Sec. IV, we perform an analysis on 7 TeV LHCb pp data, where a theory-data comparison can be found. The predictions in pPb collisions at √ s N N = 8.16 TeV, in accompanying with a discussion on the potential of determining spatial-dependent nPDFs, are presented in Sec. V. The likelihood-base approach we have used to extract σ eff,pp is documented in Appendix A. Some additional plots for the theory-data comparison at 7 TeV pp are collected in Appendix B.

II. A NEW DPS MECHANISM: FINAL STATE CORRELATIONS
In contrast to the point-like elementary particle production, there could exist another new DPS mechanism for non-point charmonium if the charm quarks are abundantly produced at short distance. As opposed to the ordinary (1), the charm and anti-charm quarks forming a J/ψ meson can come from two different partonic scatterings, where a typical Feynman diagram is sketched out at the right of Fig. 2a. In the latter configuration, the final particles J/ψ and open charm quark/hadron are strongly correlated. Such kind of channels are overlooked so far and should be in principle ubiquitous in one or multiple same flavour quarkonia production (e.g. double J/ψ [18] and triple J/ψ [19]). On the other hand, they are absent in different flavour quarkonia production (e.g. J/ψ + Υ [20]) and a quarkonium in association with other elementary particles (e.g. J/ψ + Z [21] and J/ψ + W ± [22]). We call such a contribution as DPS 2 in the following context.
By ignoring the correlation of the initial state, the expression of the differential cross section of DPS 2 in the collinear perturbative QCD factorisation and in the NRQCD approach [23] can be written as dσ DPS2 pp→J/ψ+cc = 1 σ eff,pp i,j,k,l dx 1 dx 1 dx 2 dx 2 f i/p (x 1 )f j/p (x 1 )f k/p (x 2 )f l/p (x 2 ) × n P (n) 12 dσ(ik → c( P J/ψ 2 )c(p 3 ))dσ(jl → c(p 2 )c( where the projector P (n) 12 means to cast a charm quark from one partonic scattering and a charm anti-quark from another partonic scattering into a given quantum number n = 2s+1 L [c] J , O J/ψ n is the long-distance matrix element (LDME), and the first sum runs over all possible initial partons from two beams with the corresponding parton distribution functions (PDFs) as f i/p , f j/p , f k/p , f l/p 3 . For example, for n = ß state, the precise expression is where means we have summed all the spin, colour of the external particles and have taken the average of the quantum numbers of the initial partons, andŝ ik = (p i + p k ) 2 ,ŝ jl = (p j + p l ) 2 . The phase space measure is defined as Moreover, the DPS 2 amplitude becomes where a 1 , a 2 are two colour indices, s 1 , s 2 are the (anti-)charm-quark spin components, λ is the helicity of J/ψ, is the spin projector in the nonrelativistic limit, and R J/ψ (0) is the wavefunction at the origin of J/ψ. It means in the zero correlation assumption, the DPS 2 amplitude is determined by the product of two 2 → 2 cc production amplitudes.
It is convenient to express k µ 1 ≡ 3 Without raising any ambiguity in the context, we will suppress the factorisation scale dependence in these PDFs unless when necessary.
Therefore, we have arrived at The maximally allowed phase-space integration ranges are: E2 with E 1/2 being the two beam energies, s = 4E 1 E 2 and the transverse mass of J/ψ as M T,J/ψ ≡ P 2 T,J/ψ + 4m 2 c . Due to the momentum conservation, we notice the equality p T,c = p T,c = P T ,J/ψ 2 at leading order, which is in contrast with DPS 1 , where the momentum of J/ψ is largely independent of the open (anti-)charm quark. At the LHC energies, due to the partonic luminosity, we will only consider the gluon initial state for both DPS mechanisms.

III. INITIAL CHARM CONTRIBUTION AND VARIABLE FLAVOUR NUMBER SCHEME
Another important contribution ignored so far is the initial charm contribution in SPS production. Without considering the intrinsic source of charm quarks in the proton, gg → J/ψ + cc has partially captured the initial state logarithms of the type log µ 2 h m 2 c stemming from gluon splitting into a charm quark pair as shown in the first diagram of Fig. 2b, where µ h is the typical hard scale of the process. The resummation of these logarithms to all orders in α s becomes essential when the logarithmic terms are much larger than other terms in the cross section. In such a circumstance, charm (anti-)quark should be allowed in the initial state as depicted in the middle graph of Fig. 2b (thereafter referred to as "cgF"), and the QCD DGLAP evolution of PDF automatically resum these logarithms. In the LHCb acceptance, it was indeed (roughly) estimated in Table 1 of Ref. [13] that the contribution from sea charm quark is much larger than ggF. In particular, the forward LHCb detector only measures forward particles and therefore probes the relatively high x domain of the forward beam. The initial state logarithms will be amplified with respect to the low x domain as what has been found in bottom-quark-initiated processes [24]. It is well known that the two approaches (ggF versus cgF) are compatible and complementary in phase space. The ggF calculation should be well applicable in the region close to the heavy quark threshold (i.e. µ h ∼ m c ), while the cgF calculation is better in deciphering the cross section when µ h m c . In order to correctly combine the two different calculations, a proper matching scheme must be used to subtract the overlap of the two, which is generally referring to variable flavour number scheme (VFNS). Following the pioneering work done by Aivazis, Collins, Olness and Tung (the so-called ACOT scheme) [25], the problem has been extensively studied in non-quarkonium processes in the literature [26][27][28][29][30][31][32][33][34][35][36][37][38]. Since we are only working at the lowest order for ggF and cgF processes, the matching procedure is rather simple. The only complication we have to tackle on is that since the non-relativistic nature of J/ψ, we cannot easily neglect the mass of the charm (anti-)quark in the cgF calculation. Instead, all charm quarks will be maintained as massive particles regardless where they are from. In such a situation, the differential cross section of cgF SPS c(p 1 )g( where we have taken in the initial partonic centre-of-mass system withŝ = x 1 x 2 s. 4 A similar formalism applies to gluon-charm scattering. The double counting between ggF and cgF can be attributed to first α s term of the charm PDF via with the well-known Altarelli-Parisi splitting function P qg (z) = 1 2 z 2 + (1 − z) 2 . The overlap counterterm (represented in the last diagram in Fig. 2b) that must be subtracted is Since we will be based on Pythia, the gluon splitting into a charm quark pair must be written explicitly in the event files in order to ensure the correct backward evolution of the initial state shower based on the above equation. Because of the nonzero m c , the momentum conservation and the on-shell condition cannot be guaranteed simultaneously in the splitting g → cc. A momentum reshuffling has to be imposed analogous to what has been advocated in Ref. [39] for an entirely different purpose. In the laboratory frame, the momentum of the first initial gluon is assigned as p g = x1 z E 1 , 0, 0, x1 z E 1 , and all the external momenta of cg → J/ψc have been boosted to the same frame. In particular, a boost operation has been applied to the initial charm quark p 1 → Bp 1 . Then, the momentum of anticharm is pc = p g − p 1 thanks to the momentum conservation. However, the anti-charm is not on-shell p 2 c = m 2 c without doing anything further. The momentum reshuffling is carried out as follows. We keep the invariant mass of the recoil system J/ψ + c [m 2 rec = P J/ψ + p c 2 ] invariant. The new initial gluon-gluon invariant mass square iŝ s gg = (p g + p 2 ) 2 =ŝ z = x1x2s z . In the rest frame of the new partonic system, the energies of the final anti-charm quark and the recoil system arep The three-dimensional momenta preserve the directions of the original momenta with rescalings in order to satisfy the on-shell conditions where is also guaranteed. Then, the new momenta of J/ψ and c in the same frame arẽ 4 In the centre of the hadronic collision frame P 1,2 = √ s where the boost operation B R (p) represents a Lorentz boost which turns any four-dimensional momentum to the rest frame of p. In particular, we have B R (p)p = p 2 , 0, 0, 0 . After the momentum reshuffling, all momenta of the external legs are boosted back to the lab frame.
Finally, the VFNS differential cross section is defined as where we have adopted We, therefore, arrive at a similar master formula derived in Ref. [25]. A small inconsistence may occur due to the charm quark mass. In principle, one should take the exactly same charm mass in the computations of the cross sections and in the PDF. However, we cannot keep a single m c here, because we have to stick to the pole mass in the matrix element and in the phase space, while the MS mass is usually adopted in PDF evolution. We do not bother such a small inconsistence in the paper in the view of much larger theoretical uncertainties from other sources.

A. Double parton scattering
We first compare the two DPS mechanisms (DPS 1 versus DPS 2 ) here in the LHCb 7 TeV acceptance [13]. The estimate of the DPS 1 part is different from the DPS 2 part, where the computation of the latter has been detailed in Sec. II. Since the two mesons are produced in two different partonic scatterings in DPS 1 , we opt for a data-driven approach [18,[40][41][42][43] to estimate the matrix elements of the two hard scatterings. In particular, the matrix elements of the single-inclusive prompt J/ψ and open charm hadroproduction are characterised by an empirical formula: where only the gluon-gluon initial state is retained due to its dominant partonic luminosity at the LHC, and X = g when H = J/ψ and X =C 5 when H = C. The 4 free parameters λ H , κ H , P T,H , n H are determined from the pp experimental data via a fit after convolution with a proton PDF The fitted values of these parameters for H = J/ψ, D 0 , D + , D + s and Λ + c to the LHCb double differential data d 2 σ/dP T dy [44][45][46] can be found in Table I with the CT10NLO proton PDF [47]. The factorisation scale µ F entering in the PDF is fixed dynamically as the transverse mass of the particle H. Thereby, the DPS 1 cross sections are 5C stands for the charge conjugate particle of the open (anti-)charm hadron C.
Since no correlation is considered in DPS 1 , we take an individual factorisation scale for each scattering. Finally, because LHCb [13] has presented the measured values of the cross sections for prompt J/ψ and open charm hadrons C in the acceptance, we will normalise our DPS 1 cross sections according to these measured values.
H data λH κH PT,H /GeV nH J/ψ LHCb [44,45]   Contrarily, because of the strong correlation in DPS 2 , we use the perturbative QCD approach for evaluating the scattering amplitude along with Eq.(5) (or similar variance for the colour-octet channels), where the analytic helicity amplitudes for a heavy quark pair production are known in the literature (see, e.g., in Ref. [48]). The final differential cross sections of DPS 2 are obtained via Eq. (7). The common dynamical scale µ 0 for the central value of the factorisation and renormalisation scales is H T /2, where H T is the scalar sum of the transverse mass of the final particles. The scale uncertainty for estimating the missing higher order QCD corrections is obtained by varying µ R and µ F independently around µ 0 by a factor 2. The LO proton PDF CT10LO [47] is adopted since the corresponding matrix elements are only LO accuracy here. For simplicity, we do not convolute any fragmentation function but multiply a global fragmentation fraction of It gives us the conservative upper values for the DPS 2 cross sections because the convolutions of fragmentation functions will soften the p T spectra of the open charm hadrons. The charm pole mass is taken as m c = 1.5 GeV, and the LDMEs for J/ψ and feeddown (χ c and ψ(2S)) are the Set 8 summarised in Table  3 of Ref. [49], which is originally from Refs. [50,51]. As an estimate of the order of the magnitude from DPS 2 , we restrict ourselves to consider only three cc Fock states ß,   [13]. The first quoted errors of σ DPS2 are from the scale variation and the second errors are from the PDF parameterisation.

Cross section
The comparison of the DPS 1 and DPS 2 cross sections for the 4 different final states is reported in Table II. In DPS 2 cross sections, we have quoted two errors for each cross section. The first one is stemming from the scale uncertainty and the second one is from the PDF uncertainty. The DPS 2 numbers are in anyway lower than the DPS 1 numbers by around 2 orders of magnitude. This is because in DPS 2 , a p T,C > 3 GeV kinematic cut would also imply P T,J/ψ > 6 GeV and p T,C > 3 GeV due to the strong correlation of the external momenta. In principle, the introduction of the initial primordial transverse momentum, because of the confinement and the uncertainty relation, can violate such a relation. However, we have explicitly checked that the kicks of the initial state in the transverse plane by a few GeV can maximally enhance the DPS 2 cross sections by a factor of ∼ 5, which are still at least one order of the magnitude smaller than DPS 1 . Such a conclusion will not change if we take other LDME sets in Ref. [49]. It is, however, not necessary to be still true that DPS 2 will be negligible compared to DPS 1 under other conditions. We will examine it again in proton-lead collisions at √ s N N = 8.16 TeV in Sec. V, where there is no lower p T cut imposed.
B. Single parton scattering

Colour-octet and feeddown contributions
A brief survey of the possible impact of the colour-octet contribution [52] and of χ c feeddown [2] reveals that they only alter the dominant colour-singlet ß predictions marginally until reaching to the tail of the P T,J/ψ spectra (e.g. P T,J/ψ > 15 GeV), which are outside the LHCb kinematic requirement. Here, we will assess their contributions with the diverse LDME fits on the market by considering ggF SPS. The simulation is carried out by the joint usage of HELAC-Onia 2.0 [14,15] and Pythia 8.186 [16] with CT10LO as the proton PDF. In order to match what has been implemented in Pythia 8, the charm quark mass m c is taken the half of the physical mass of the meson in the colour-singlet channels (i.e., ß for J/ψ and ψ(2S), and J for χ cJ ), while a 100 MeV increment is applied to m c for the other S-and P-wave colour-octet channels (i.e., J for J/ψ and ψ(2S), and 3 S [8] 1 for χ cJ ). Totally, we have considered 18 Fock state channels. Each channel consists both gluon-gluon and quark-antiquark initial states. The central scale is same as DPS 2 case, i.e., H T /2. The values of LDMEs from worldwide data fits are extremely disparate and inconsistent among different groups. The status has been well summarised in Sec. 5.2 of Ref. [49]. We take a global survey in order to minimise the possible bias introduced by a concrete LDME fit. A summary of 9 LDME sets is reported in Table 3 of Ref. [49], which are originally taken from Refs. [50,51,[53][54][55][56][57][58][59]. We follow the same notation here and refrain from tabulating them again. We have reported the ratios of the prompt cross sections over the colour-singlet direct production J/ψ 6 (i.e. no colour-octet transition and no feeddown) in Table. III. For a given LDME set, the ratios only mildly depend on the species of the open charm hadron. On the other hand, the dependencies of such ratios on LDMEs are more significant. Taking J/ψ + D 0 as an example, the ratio ranges from 1.06 (Set 4) to 1.65 (Set 2). This corroborates the similar conclusion drawn based on earlier (partial) studies: colour-octet and feeddown contributions are not able to dramatically enhance the cross section.

Fixed flavour number scheme versus variable flavour number scheme
We compare the fixed flavour number scheme results and the variable flavour number scheme results in SPS. Again, for the sake of simplicity, we only consider the direct J/ψ production via the leading ß channel since the inclusion of the colour-octet and feeddown channels has been discussed in the previous subsection. Similarly to the setup in the last subsection, the LDME O J/ψ ß has been fixed as 1.16 GeV 3 , and the central scale µ 0 = H T /2. The simulations are carried out by jointly using HELAC-Onia 2.0 [14,15] and Pythia 8.235 [60] with primordial k T and underlying event enabled. At variance with what has been done before, we consider several proton PDF choices here. There are one leading order (LO) PDF (CT10LO [47]), 4 next-to-leading order (NLO) PDF sets (CT10NLO [47], CT14NLO [61], MMHT14NLO [62] and NNPDF3.1NLO [63]) and a next-to-next-to-leading order (NNLO) PDF set (CT14NNLO [61]). Although all of them are global-fitted PDFs based on variable flavour number schemes, there are a few significant differences among them. Notably, they differ by: 1) perturbative orders in the scale evolution of PDF and α s and in matrix elements to fit the PDF; 2) PDF parameterisations, statistical fit methodologies and uncertainty estimations; 3) values of parameters (e.g. α s (m 2 Z ) and heavy quark masses) and input experimental data; 4) the concrete implementation of variable flavour number schemes. None of them has been introduced the non-perturbative source of charm and bottom (anti-)quark densities.   (18)] separately. Two theoretical errors are quoted for each configuration along with its central value. The first one is from the standard 9-point renormalisation and factorisation scale variation µ R/F = ξ R/F µ 0 , ξ R/F ∈ {1.0, 0.5, 2.0}. The second error is estimated from 68% confidence level (CL) PDF parameterisation uncertainty. In general, because we are working at the lowest order (i.e. LO) in the strong coupling α s for the matrix elements and at rather low scales (a few GeV), the scale uncertainties anyway dominate the theoretical uncertainties. They shift the predicted cross sections up and down by factors of more than 3 around the central values. The effect is more striking in processes with only three charm (anti-)quarks involved in the hard scatterings (i.e. cgF, CT) , which is anticipated from their even lower scales than ggF, where in the latter there are four heavy (anti-)quarks. Such theoretical errors can only be systematically reduced by including higher-order QCD terms in the α s perturbative series. Charm quark hadroproduction up to NNLO reveals that the precise predictions with higher-order α s corrections usually lie at the upper limits of the corresponding LO predictions. Such a statement can of course be altered by considering a different process and by choosing a different central scale. Given the absence of higher-order QCD calculations of the J/ψ plus open charm process nowadays, we allow the missing higher-order corrections can be any values allowed by the scale variation, and the differential shapes can be distorted arbitrarily within the scale uncertainty bands. A remarkable observation from Table IV is that the original ggF calculations in 3 flavour number scheme (3FS) proposed in the literature [1,4,5] largely underestimate the SPS predictions of VFNS. VFNS enhances 3FS ggF central predictions by a factor 3, while the upper limits due to the scale variation are more than 4 times larger after including the (anti-)charm initial state. This implies the importance of resumming the initial collinear logarithms α n s log k µ 2 F m 2 c with 0 < k ≤ n in our interested kinematic regime. Hence, the statement of "SPS can be negligible" should be revised, and it sheds light on resolving the tensions between the LHCb data and the theoretical predictions discussed in the introduction. On the other hand, the PDF uncertainty is subdominant compared to the scale uncertainty. It introduces additional ∼ 15% and ∼ 30% theoretical errors for ggF and VFNS results. It is, however, interesting to notice that not all PDF set gives reasonable estimates of the true PDF uncertainty.
C. σ eff,pp and theory-data comparison

Determination of σ eff,pp
Following the prescription detailed in Appendix A, we redetermine σ eff,pp from the LHCb data via the likelihoodbased approach. In order to maximise our prediction power and to minimise the possible bias, we only choose the integrated cross section, the normalised shapes of P T,J/ψ and invariant mass M (J/ψ + D 0 ) distributions from the LHCb measurement of J/ψ +D 0 final state. We have checked that the integrated cross section alone does not give any meaningful constraint for the effective cross section σ eff,pp , which is plagued with huge SPS theoretical uncertainties (especially the scale uncertainty) as shown in the Sec. IV B. Not all data in the aforementioned distributions will be adopted in the likelihood fit. We only select the 6 P T,J/ψ ∈ [1.0, 4.0] GeV data and 10 M (J/ψ + D 0 ) ∈ [5.5, 10.5] GeV data. The very low P T,J/ψ data (P T,J/ψ < 1 GeV) are excluded because of the imperfection in our DPS modelling in the regime (cf. Fig.1a and 1b in Ref. [43]). For the same reason, the first near-threshold bin 5.0 < M (J/ψ+D 0 )/GeV < 5.5 has been excluded as well. The reasons for not taking into account the data in the tails of the two spectra are twofold. The statistical Monte Carlo fluctuations in the theoretical curves start to be significant. In the P T,J/ψ distribution, we also encounter an unexpected large fluctuation in bin 4.0 < P T,J/ψ /GeV < 4.5 of the LHCb data, which is beyond the size of the given experimental error. In total, after accounting for the integrated cross section, we have used 17 experimental data in our likelihood fit procedure to extract the value of σ eff,pp .
In Fig. 3, we show the marginalised constraints on σ eff,pp from the 7 TeV LHCb J/ψ + D 0 data. The posterior probability density is defined as L 17 (σ eff,pp )  Fig. 4. For VFNS SPS with various PDF sets as reported in Table IV, the best-fitted values of σ eff,pp span from 29.9 mb to 34.8 mb. The 68% quantiles of the effective cross section are reported both in the figure legend of Fig. 3 and in Fig. 4. The worst precision is from the SPS with the NNPDF3.1NLO PDF, but is still better than 12%. All of the inferred values of σ eff,pp from various PDF sets are consistent with each other, which are within one standard deviations in Fig. 4. The inclusion of VFNS SPS dramatically increases the σ eff,pp values by a factor 2 with respect to those assuming negligible SPS in Ref. [13]. As we will show in the later subsections, the theoretical calculation and the LHCb data will also be largely reconciled. Given the PDF dependence is marginal with respect to the other sources of theoretical errors, in the following, we will only use CT14NLO in the SPS results. The final determination of σ eff,pp with this PDF is 34.8 +1.2 −2.5 mb, posing the most precise (around 7%) constraint. Such a σ eff,pp value from J/ψ + D 0 will also be applied to other three final states J/ψ + C(C = D + , D + s , Λ + c ), because they basically undergo the same hard scattering process 7 . With this setup, we start to have our theory-data comparison.

Integrated cross sections
After combining all the aforementioned information, we are able to refine the theoretical calculations, in particular the SPS postdictions. The comparison to the LHCb measurement can be found both in Table V and in Fig. 5. We call the coverage of the prompt over colour-singlet J/ψ ratios listed in Table III as "LDME uncertainty". The central values are rescaled with the ratio numbers of the LDME Set 5, which should be taken with a grain of salt. For the sake of being conservative, we take the envelope of the 9-point scale variation in 9 LDME sets as our combined scale+LDME uncertainty estimate. PDF uncertainty will be either given separately (e.g. in Table V)  way with other theoretical uncertainties. As pointed out in the last subsection, we use σ eff,pp = 34.8 mb in the following context in DPS estimates. Since our DPS results are more-or-less data driven (modulo σ eff,pp ), we do not associate any theoretical uncertainty to them. As announced, the scale uncertainty is the dominant theoretical error, which can be clearly seen in Fig. 5. A higher-order QCD calculation, at least NLO, is demanding in the near future to improve the accuracy. Nevertheless, assuming that the scale variation has successfully captured the correct size of the missing higher-order α s terms, the VFNS SPS can be as sizeable as the experimental data. Therefore, from the viewpoint of the integrated cross section, in principle, no DPS is needed. However, the current prominent uncertainty cannot forbid a significant DPS contribution. We will also give a remark on the open charm baryon Λ + c here. The baryon Λ + c has been shown at the LHC that the traditional single-parton fragmentation function is insufficient to describe the single-inclusive data [64], which calls for novel mechanisms (e.g. the coalescence mechanism [65]). Hence, one should be careful to interpret the J/ψ + Λ + c data here, since we only include the fragmentation contribution in our SPS simulations.   [13]. The quoted theoretical errors in the second and fourth columns represent scale+LDME and PDF uncertainties, respectively. The first (second) uncertainty in the LHCb data column is statistical (systematic).

Differential distributions
We now turn to the differential distributions, which are originally in tensions with the DPS-dominance hypothesis. Given the similarity among different open charm hadron species, we will only focus on J/ψ + D 0 here, while the comparison for other hadron species can be found in Appendix B. Following what has been given in the LHCb paper [13], the comparison is only performed at the normalised distribution (shape) level, i.e. the sum of all bins gives unity. This will help to cancel global and correlated systematical experimental errors, which are present in the absolute differential distributions. On the contrary, we view the scale uncertainties in SPS are not bin-by-bin correlated. In fact, the missing higher order quantum radiative corrections may distort the shapes within the given uncertainty bands. For other theoretical uncertainties, they can be treated in the full correlation way. However, they are subdominant and unimportant from this perspective.
The theory and data comparison for the five distributions is presented in Fig. 6. The distributions are the transverse momenta of J/ψ (Fig. 6a) and D 0 (Fig. 6c) mesons, the invariant mass of the two mesons (Fig. 6b), the rapidity gap ∆y(J/ψ, D 0 ) = y J/ψ − y D 0 (Fig. 6d) and the azimuthal angle difference ∆φ(J/ψ, D 0 ) = φ J/ψ − φ D 0 (Fig. 6e) between the meson pair. We have marked the data those enter into our likelihood fit in determining σ eff,pp in green. In general, the previous tensions are greatly alleviated. The reasonable agreements have been achieved in most of the distributions, except the first two bins in the P T,J/ψ distribution. The disagreement in the first two bins is anticipated due to the limitation of our DPS modelling. Although there are a few random fluctuations in the experimental data, almost all of them lie within the grey bands.
The spectra of P T,J/ψ and M (J/ψ + D 0 ) are quite distinct between SPS and DPS. This fact reflects that J/ψ undergoes the single-charm fragmentation process in J/ψ + c production at high P T , while several competitive mechanisms are present in DPS J/ψ production, which is equivalent to the single inclusive J/ψ process. The LHCb data favour the hard spectra close to the SPS ones, implying the necessities of the SPS events in interpreting the data. On the other hand, SPS curves are quite similar to the DPS ones in the p T,D 0 and ∆y(J/ψ, D 0 ) distributions. In the p T,D 0 distribution, D 0 are generated via the charm quark fragmentation in both SPS and DPS. Its shape, hence, has no discrimination power to separate the two contributions. In the rapidity gap ∆y(J/ψ, D 0 ) distribution, DPS shape turns out to be very symmetric around ∆y(J/ψ, D 0 ) = 0 because the two mesons are mainly produced in two independent scattering subprocesses, while due to the correlation in SPS, the D 0 meson is slightly more forward than the J/ψ meson. However, such an asymmetric behaviour in the SPS distribution is still quite mild. The azimuthal angle difference ∆φ(J/ψ, D 0 ), in principle, could be useful to disentangle the SPS and the DPS mechanisms, because DPS contribution is flat due to the dominance of DPS 1 , while the SPS has no reason to be flat due to the correlation of the final particles. However, as pointed out several times before, the nonflat ∆φ(J/ψ, D 0 ) distribution can be diluted by the introduction of the initial k T smearing. The smearing effect can make the original nontrivial structures in ∆φ(J/ψ, D 0 ) fading out. In our simulation, such an effect is modelled by both the initial state radiation and the primordial k T in the beam remnants implemented in Pythia8. However, we should remind the reader that the primordial k T is pure phenomenological and suffers from the uncertainty in model dependence. With our simulation, it seems that the shape of the SPS is inadequate to describe the experimental data of ∆φ(J/ψ, D 0 ), while the corresponding DPS shape matches the data in a better way. The SPS contribution peaks at the away side ∆φ(J/ψ, D 0 ) = π, i.e. back-to-back in the transverse plane. The sum of the SPS and the DPS has a right shape, which indicates that the DPS process should not be ignored as well.

V. PREDICTIONS FOR PROTON-LEAD COLLISIONS AT
√ sNN = 8.16 TEV After understanding the pp data, we are in the position to present our theoretical predictions for the upcoming LHCb measurement in proton-lead (pPb) collisions at the average nucleon-nucleon center-of-mass energy √ s N N = 8.16 TeV. 8 As advocated in Ref. [11], such a measurement will be very useful, because it will, for the first time, reveal both the nucleus geometrical enhancement effect [67] and the impact-parameter dependence of the nuclear modification on the incoming parton flux. In this context, we assume the only dominant cold nuclear matter effect has been encoded in the universal nuclear PDFs (nPDFs). Such a hypothesis has been effectively justified for (non-excited) heavy flavour production at the LHC [12,43], where the nPDF effect is quite significant and therefore must be taken into account. In order to quantify the spatial dependence of nPDFs, we assume such a dependence is only related to the local nuclear thickness function T A ( b ), where b is the two-dimensional impact parameter vector in the colliding transverse plane. The nuclear modification R A k (x, b ) for a parton k = g, q,q at the position b is expressed as where R A k (x) is simply the ratio of nPDF over free-nucleon PDF for the parton k at the Bjorken x. G() could be any function that satisfies the condition We choose a test function a , and opt for the simple hard-sphere form for the thickness function T A . 9 After incorporating both the nuclear-collision geometry and the spatial-dependent nPDFs, we have the following DPS i differential cross section in proton-nucleus collisions [11]: dσ DPSi pA→J/ψ+C = Adσ DPSi J/ψ+C,11 where dσ DPSi J/ψ+C,kl is the DPS i (differential) cross section in (free or bounded) nucleon-nucleon collisions. The convoluted PDFs in dσ DPSi J/ψ+C,kl , for the two initial partons from the ion beam, are two (spatial-averaged) nPDFs when kl = 11, one nPDF and one free-nucleon PDF when kl = 10 or 01, and two free-nucleon PDFs when kl = 00, respectively. It means dσ DPSi J/ψ+C,00 = dσ DPSi pp→J/ψ+C without considering the possible isospin effect (e.g. the gluon initial state in our case). In this context, we take EPPS16 [68] with the LHC J/ψ or D 0 constraint derived in Ref. [12] as our nPDF set, 10 and CT14NLO as the corresponding proton PDF. According to our pp fit, σ eff,pp is fixed to 34.8 mb, and A = 208 and R A = 6.624 fm for Pb. In the remainder of the paper, without losing generality, all the kinematics will be defined in the centre-of-mass frame of the averaged per nucleon-nucleon collision.

A. Cross sections
We have reported the breakdown of cross sections from three different sources, VFNS SPS, DPS 1 and DPS 2 , within the LHCb detector acceptance [69] in Table VI, where the power a appearing in the G() is fixed to the widespread value a = 1, in which the nuclear modification is proportional to the thickness function. The forward (backward) region imposes both the J/ψ meson and the open (anti-)charm hadron C being in the rapidity interval 1.7 < y J/ψ , y C < 3.7 (−4.7 < y J/ψ , y C < −2.7). 11 Several theoretical uncertainties are associated to various contributions. In the catagory of VFNS SPS, the first quoted errors are from the 9-point scale variation and the LDME uncertainty, while the second errors stand for both the nPDF and proton PDF parameterisation uncertainties. The two errors in DPS 2 numbers are the scale uncertainty and the reweighted nPDF uncertainty, where, as mentioned in the subsection IV A, we have only included three Fock states (ß,  Table 3 in Ref. [49]. Moreover, only 68% nPDF uncertainties are associated to the DPS 1 cross sections. Similar as the finding in the 7 TeV pp case, DPS 2 contributions are at least 20 times smaller than DPS 1 even taking into account the large scale uncertainty. On the other hand, the upper limits of the VFNS SPS predictions are comparable with the DPS 1 components. Since there is no lower p T cut, we are probing even lower scales than the 7 TeV pp case. The scale uncertainty in VFNS SPS is extremely large. Such a situation is inevitable before carrying out a computation including higher-order QCD corrections. For the most abundant process J/ψ + D 0 , the pPb cross sections are predicted to be around 0. a that quantifies the impact-parameter dependence of the nuclear modification via Eq. (25). There are predictions both in the forward (left) and in the backward (right) rapidity intervals of the LHCb detector.
in the both fiducial regions. Given the similarities of the results in various open charm hadron species, we will only concentrate on J/ψ + D 0 plus its charge conjugate mode in the following. The potential in constraining impact-parameter dependent gluon nPDF from the process is very intriguing. The cross sections of J/ψ + D 0 are displayed in the functions of a ∈ [0.0, 4.0] in Fig. 7. As obvious, the VFNS SPS (red bands) are independent of the a value. 12 The a dependencies of the two DPS contributions (DPS 1 versus DPS 2 ) are more significant in the forward region than in the backward region. This is simply because the nuclear modification in the latter is much smaller. Fig. 7a tells us if the impact-parameter dependence of the gluon nuclear modification is strong enough (a > 3), we can easily observe the abundant J/ψ + D 0 events in the forward detection. In turn, this parameter space can be excluded if the measured cross section is smaller than 1 mb. The challenging case is of course ) is weaker, in other words when a ≤ 3 in our simplest parameterisation of G(). Aside from the integrated cross sections, the differential distributions can provide additional useful information. The invariant mass M distributions of the J/ψ + D 0 system can be found in Fig. 8 Both forward (left) and backward (right) regions are considered.
forward (backward) rapidity region. The two DPS components have been summed, in which DPS 1 always dominants. As our showcases, we consider five different "a" values in our G() parameterisation. They are a = 0 (zero spatial dependence, green filled circle), 1 (most widely used assumption, blue empty circle), 2 (purple empty square), 3 (brown filled square) and 4 (orange). The same layout will be applied for the other distributions too. In order to make the plots less busy, only one DPS uncertainty band a = 1 is associated to the DPS predictions in each plot. In the forward region, the probing x in nPDF resides in the shadowing region, which has the significant nuclear modification effects. On the contrary, the backward region is exploring the x regime close to the transition from the shadowing to the antishadowing. Hence, the nuclear modification is smaller. Likewise the observation in the total cross sections, the DPS a ≤ 3 distributions are quite similar and close to each other, while the very different features can be observed when a = 4. In the forward case (Fig. 8a), the a = 4 DPS distribution is strongly enhanced uniformly in all bins with respect to the a = 3 distribution. On the other hand, the backward a = 4 curve in Fig. 8b features a more peculiar shape. In particular, the DPS prediction becomes negative if 5.5 < M (J/ψ + D 0 )/GeV < 8.5. This can be understood from Eq. (27) The differential cross section would be negative when the second piece in the brackets is larger than the other two terms. Such a pathological behaviour, however, could be very sensitive to the fine tuning among the different input parameters. Moreover, although the SPS predictions are plagued with sizeable theoretical uncertainties, the DPS cross section, whatever the value of a, starts to overshoot the VFNS SPS cross section in large M (J/ψ + D 0 ) in the backward region, especially when the mass is from 10 GeV to 18 GeV. It benefits from the reduced scale uncertainty in the VFNS SPS cross section at large scale. Such a specific kinematic region provides a useful mean to improve the purity of the DPS signal. Another kind of interesting distributions is the transverse momentum P T of the two meson pair, as shown in Fig. 9. Due to the importance of initial (anti-)charm, the P T (J/ψ + D 0 ) spectra of VFNS SPS are soft, because the produced J/ψ meson and the (anti-)charm quark are back-to-back in the LO 4FS calculation. Such a picture would not change even with higher order QCD corrections, since the two mesons are still likely along the beam pipe axis following the fragmentation of the (anti-)charm quark in the new topology appearing in the real emission diagrams. As opposite, the two mesons in the DPS events are almost independently produced. Hence, the DPS curves in Fig. 9 are generally harder than the VFNS SPS curves. A lower P T (J/ψ+D 0 ) cut can substantially enhance the DPS fraction in the events in the both interested rapidity intervals. P T (J/ψ + D 0 ) could be used to define the control regions of the DPS events on the experimental side. Similar to the invariant mass distributions, the strongest impact-parameter hypothesis with a = 4 yields a peculiar dip structure in Fig. 9b, which again can be attributed to the delicate cancellation among the The rapidity gap ∆y between the two detected mesons is widely used to differentiate the DPS and SPS contributions in the double J/ψ production process. Here, we examine the same quantity for the J/ψ + D 0 associated production process. The corresponding distributions at the two rapidities can be found in Fig. 10. In order to fit in the frame, we have multiplied a factor 0.05 to the DPS a = 4 curve in Fig. 10a. because the two tagged mesons are not identical, these distributions are not necessary to be symmetric around ∆y(J/ψ, D 0 ) = 0. This has indeed been observed in the VFNS SPS distributions. In both Figs. 10a and 10b, our simulation shows that J/ψ tends to be more central in rapidity than the D 0 meson. Due to the zero correlation of the two mesons in the DPS 1 events, the DPS distributions are (very close to be) symmetric around ∆y(J/ψ, D 0 ) = 0. A small asymmetric feature in Fig. 10b simply reflects some delicate details of our DPS 1 template, which is driven by the single-inclusive particle production data in pp collisions. This asymmetry is amplified further in the DPS a = 4 curve. The aforementioned cancellation results in a very narrow ∆y(J/ψ, D 0 ) distribution when a = 4 and the two mesons are backward.
The last distributions, which we are interested in, are the azimuthal angle difference ∆φ between J/ψ and D 0 . They are exhibited in Fig. 11. The ∆φ(J/ψ, D 0 ) observable is known to be very sensitive to the transverse momentum of the incoming partons from the intrinsic nonperturbative source and the perturbative initial shower. These effects have been taken into account in the VFNS SPS simulation via Pythia8 and in the template of DPS 1 , while they are absent in the LO calculation of DPS 2 . The peak structures at ∆φ(J/ψ, D 0 ) = π in Fig. 11 are from the DPS 2 events. We have examined that such peaks quickly fade out if we allow nonzero k T of the initial partons. After taking into account such an initial-k T smearing effect, the final DPS events are uniformly distributed in all bins, while the VFNS SPS events tend to populate at the back-to-back region ∆φ(J/ψ, D 0 ) π. However, the distinction between the two is still not remarkable enough given the current theory precision.

B. Nuclear modification factors
A traditionally useful observable in proton-nucleus collisions is the nuclear modification factor, which quantifies the relative modification of the cross section in pA collisions with respect to its counterpart in pp collisions. For our interested process, it is defined as We have used the abbreviation r eff,A ≡ (A − 1) . For lead Pb, r eff,A 5.23 σ eff,pp 34.8 mb . It means that the terms linear in r eff,A are more important in the nuclear modification factor R DPSi pA→J/ψ+C unless σ eff,pp is smaller than 7 mb. Thus, the DPS R pPb should be more sensitive to the value of σ eff,pp than the absolute cross sections in pPb collisions. In the following context, we will also consider the uncertainty in σ eff,pp = 34.8 +1.2 −2.5 mb in evaluating the DPS R pPb values. We want to point out that the nuclear modification factor R pPb also provides us a new handle to simultaneously extract σ eff,pp and G(). A caveat is that since the fiducial volumes of the 7 TeV pp and 8.16 TeV pPb data are not completely identical, in principle, the identify of σ eff,pp in the two measurements may not be fully justified if a slight kinematic dependence of σ eff,pp is allowed for instance. Therefore, a dedicated pp measurement under the same or at least similar condition would help to reduce such an uncertainty.
The R pPb values in the two LHCb fiducial regions are reported in Table VII. We only consider J/ψ + D 0 in this section since other three final states share the very similar nuclear modification factors. R pPb predictions from VFNS SPS are shown in the second row, while those from DPS, as well as its breakdown into DPS 1 and DPS 2 parts, can be found in the last, third and fourth rows respectively. The R pPb of VFNS SPS is smaller than unity in the forward rapidity interval 1.7 < y J/ψ , y D 0 < 3.7. Such a fact is anticipated because we are probing the shadowing region of nPDF. Alternatively, the parton nuclear modification encoded in nPDF is very small when −4.7 < y J/ψ , y D 0 < −2.7, which yields the SPS R pPb very close to unity. Several theoretical uncertainties are quoted. The first errors, including those in the parenthesis, are from the scale variation. The usual 7-and 9-point (in the parenthesis) scale variation 14 errors are very large when 1.7 < y J/ψ , y D 0 < 3.7. The big scale uncertainty stems from the low factorisation scale µ F = 0.5µ 0 in R Pb k (x, µ 2 F ) with k = g, c,c, which could be close to the charm-quark mass threshold. In the later analysis of the differential R pPb , we will explicitly show that such a scale dependence actually rapidly diminishes when the hard scale increases. In the backward region, because R Pb k (x, µ 2 F ) ∼ 1, the scale uncertainty largely cancels in the ratio and it becomes marginal. The second quoted errors are from the R Pb k (x) parameterisation at the initial scale. They have been greatly constrained by the inclusive heavy-flavour data in Ref. [12]. We have taken the envelope of the three reweighted EPPS16 grids, which have utilised the single heavy-flavour meson production data in three different factorisation scale choices. In addition, we have assumed that there is no correlation of the hard scales between the SPS J/ψ + D 0 production and the single inclusive meson production processes. The final nPDF uncertainty is 18% for 1.7 < y J/ψ , y D 0 < 3.7 and 6% for −4.7 < y J/ψ , y D 0 < −2.7.   We have only shown the DPS R pPb 's with a = 1 in Table VII, while their a dependencies in a wider range a ∈ [0, 4] can be found in Fig. 12. Due to the geometrical enhancement effect from those terms proportional to r eff,A , the DPS R pPb 's are generically much larger than unity whatever the value of a, as long as σ eff,pp is not too small. With our nominal value of σ eff,pp from the pp data, we have R DPS pPb ∼ 3 for the forward mesons and ∼ 6 for the backward mesons when a ≤ 3. A factor 2 suppression between the two rapidity intervals is attributed to the double nuclear parton suppression R Pb k in R DPS J/ψ+D 0 ,11 . In the breakdown rows (DPS 1 and DPS 2 ) of Table VII, only the σ eff,pp (first errors) and nPDF (second errors) uncertainties are considered. Other commonly present errors both in pPb and in pp can be taken as zero because of the cancellations in the ratio. As opposite, the cancellation is not guaranteed when we sum the two components of DPS concerning the scale variation in DPS 2 . Therefore, the last row of Table VII shows three different errors. The first ones are because of the 9-point scale variation in DPS 2 , and the second (third) ones are from the uncertainty in σ eff,pp (in nPDF). Unlike VFNS SPS, the dominant uncertainty in DPS is due to the nPDF parameterisation. Such an uncertainty changes the central value by 15% (forward) and 8% (backward). When a > 3, R DPS pPb increases dramatically in the forward region (see Fig. 12a). It, therefore, provides a smoking gun signal, if one is able to observe an anomalously large nuclear modification factor with an experiment. Finally, we comment on an additional feature in the layout of the figures in this subsection. There are two bands associated to each VFNS prediction: one dark-red band and one light-red band. The dark (light) bands are the combined theoretical uncertainty with the 9-point (7-point) scale variation. The differences between the two will be quite significant at low scales but will not be even visible if we go to higher scales.
We now turn to the discussions of the kinematical distributions of R pPb . Following the absolute cross section discussions, we consider 4 different observables: the invariant mass M (J/ψ + D 0 ) (Fig. 13), the transverse momentum P T (J/ψ + D 0 ) (Fig. 14), the rapidity (Fig. 15) and azimuthal angle (Fig. 16) From Fig. 15b, if the J/ψ and D 0 mesons are more asymmetrically produced in rapidity, the discrimination power for the different curves in a of the DPS events becomes larger.
Since our predictions are plagued with the very large theoretical uncertainties, in particular in SPS, we do not show the combined R pPb after encompassing both SPS and DPS contributions. Before closing this section, we would like to comment on how to obtain R pPb based on our results. There are at least two ways. If the fraction f pPb DPS of DPS events in pPb collisions can be better known (e.g. in some control regions), we can obtain the combined R pPb via On the other hand, if one is able to determine the DPS fraction f pp DPS in a corresponding pp measurement, ideally under the same condition, we are able to derive R pPb via VI. CONCLUSIONS In this paper, we have carried out a thorough study on the process of a J/ψ meson production in association with an open (anti-)charm hadron at the LHC. As our first step, we have tried to understand the existing LHCb measurement in pp collisions at √ s = 7 TeV, where there were apparent tensions between the DPS dominance picture and the experimental data. We proposed a new DPS mechanism and revised the SPS Monte Carlo simulation based on HELAC-Onia 2.0+Pythia8 in order to resolve the discrepancies. The novel DPS mechanism DPS 2 has peculiar features with respect to the conventional DPS procedure DPS 1 . The former is strongly correlated among the particles in the final state, while the latter is expected to correlate only in a weak way between the two distinct scatterings. In our specific process, DPS 2 , however, turns out to be small. Nevertheless, such a new DPS 2 process is ubiquitous for processes involving more-than-one pairs of same-flavour heavy quarks with at least a quarkonium, like J/ψ + cc and J/ψ + J/ψ etc. The DPS 2 contributions in the processes other than J/ψ + cc deserve being investigated in the future. On the other hand, in order to scrutinize the LHCb data, we have performed a proper matching between 3FS and 4FS calculations in the SPS simulation. Our calculations demonstrate that the 4FS part, which resummes the large logarithms of initial gluon-to-charm splitting, is indispensable to account for the robust SPS predictions and also the LHCb data. The matched VFNS calculation significantly enhances the SPS cross sections and alters the DPS-dominance conclusion, which was based on the 3FS-alone simulation. To the best our knowledge, it is also a first VFNS calculation for a quarkonium process in the literature. After encompassing all components, we (re)determine the effective cross section σ eff,pp entering into the DPS formula from the LHCb data with a likelihood-based approach. Without big surprise, the new determinations yield σ eff,pp 30 mb with 6 different proton PDFs, which are a factor 2 larger than the previous fit based on the DPS-alone hypothesis. These improvements allow us to fill the gaps between theory and experiment, albeit with large theoretical uncertainties. This clearly calls for an improved simulation in the future that contains higher-order quantum corrections in order to ensure whether the discrepancies have really been gone.
Furthermore, we have also presented our predictions for the same process but in pPb minimum-bias collisions at √ s N N = 8.16 TeV. The corresponding experimental measurement is carrying on by the LHCb collaboration. We identify a few kinematical regions that would be very useful in improving the purity of the DPS events. Namely, they are the large M (J/ψ + D 0 ) and high P T (J/ψ + D 0 ) regions. We also exploited the potential to extract the spatialdependent nPDF information from this process. The nuclear modification factors R pPb are very different between DPS and SPS, and are also significantly depending on the (unknown) impact-parameter dependence of nPDF. Therefore, with enough statistics, the LHC can definitely tell us something on how initial partons are impacted when they are located at different positions in a nucleus. We are looking forward to see the emergence of more measurements at the LHC in the near future.
In particular, the central value of σ eff,pp is determined by maximising L (σ eff,pp ), which we denote as σ eff,pp = δ 0 , and its 68% CL error is determined via The final 68% confidence interval is σ eff,pp = δ 0 +δ+ −δ− . In general, δ + and δ − should be not identical, and σ eff,pp poses an asymmetric error.
In the presence of several (differential) cross section measurements, we have to opt for an approach to combine the individual determinations of σ eff,pp . We adopt a likelihood-based approach as advocated in Refs. [71,72], in which a global likelihood function is constructed from the probability distribution functions of individual extractions. For simplicity, we do not breakdown the errors of individual determinations into different sources and take all of them uncorrelated. Such a simple and conservative treatment allows us to avoid introducing additional nuisance parameters.
Let us assume that we want to combine n different determinations with the ith one as δ where The global likelihood function is then and the test statistic is q n (σ eff,pp ) = −2 log L n (σ eff,pp ) L n (σ eff,pp ) , where L n (σ eff,pp ) is maximised when σ eff,pp =σ eff,pp . The test statistic q n is always positive or zero. The 1σ (i.e. 68%) confidence interval is extracted from q n = 1.