Double D-meson production in proton-proton and proton-lead collisions at the LHC

TL;DR

This paper models the simultaneous production of two D mesons in proton-proton and proton-lead collisions at the LHC, using next-to-leading order QCD, and predicts measurable double-parton scattering signals in proton-lead collisions.

Contribution

It introduces a QCD-based framework for double D-meson production that incorporates double-parton scatterings and makes predictions for proton-lead collisions at the LHC.

Findings

Framework agrees with LHCb proton-proton data

Predicts measurable double D-meson production in proton-lead collisions

Provides evidence for double-parton scattering at perturbative scales

Abstract

We consider the simultaneous production of two heavy-flavoured hadrons - particularly D mesons - at the LHC. We base our calculations on collinearly factorized QCD at next-to-leading order, using the contemporary parton distribution functions and D-meson fragmentation functions. The contributions of double-parton scatterings are included in the approximation of independent partonic interactions. Our framework benchmarks well with the available proton-proton data from the LHCb collaboration giving us confidence to make predictions for proton-lead collisions. Our results indicate that the double D-meson production in proton-lead collisions should be measurable at the LHCb kinematics with the already collected Run-II data, and should provide evidence for double-parton scattering at perturbative scales with a nuclear target.

Figures

Equations

$$\frac{\mathrm{d}{\sigma }_{AB\to a+b+X}}{{\mathrm{d}}^3{\vec{p}}^a{\mathrm{d}}^3{\vec{p}}^b}=AB\left[\frac{\mathrm{d}{\sigma }_{nn\to a+b+X}^{\mathrm{sps}}}{{\mathrm{d}}^3{\vec{p}}^a{\mathrm{d}}^3{\vec{p}}^b}+\frac{m}{{\sigma }_{\mathrm{eff}}^{AB}}\frac{\mathrm{d}{\sigma }_{nn\to a+X}^{\mathrm{sps}}}{{\mathrm{d}}^3{\vec{p}}^a}\frac{\mathrm{d}{\sigma }_{nn\to b+X}^{\mathrm{sps}}}{{\mathrm{d}}^3{\vec{p}}^b}\right]$$

$$\displaystyle\frac{1}{\sigma_{\rm eff}^{AB}}\equiv\Bigg{\{}\frac{1}{\sigma_{\rm eff}}$$

$$+\frac{(A-1)}{A^2}\int {\mathrm{d}}^2\vec{B}{\left[T_{nA}\left(\vec{B}\right)\right]}^2$$

$$\displaystyle+\frac{(A-1)(B-1)}{(AB)^{2}}\int\mathrm{d}^{2}\vec{B}\left[T_{AB}\left(\vec{B}\right)\right]^{2}\Bigg{\}}\,.$$

$$\frac{1}{{\sigma }_{\mathrm{eff}}}=\int {\mathrm{d}}^2\vec{b}{\left[t_{\mathrm{nn}}\left(\vec{b}\right)\right]}^2,$$

$$t_{\mathrm{nn}}(\vec{b})\equiv {\int }_{-\infty }^{\infty }{\mathrm{d}}^2\vec{s}{t}_{\mathrm{n}}(\vec{s}+\vec{b}/2)t_{\mathrm{n}}(\vec{s}-\vec{b}/2),$$

$$t_{\mathrm{n}}(\vec{s})\equiv {\int }_{-\infty }^{\infty }\mathrm{d}z{\rho }^{\mathrm{n}}(\vec{s},z).$$

$$T_{\mathrm{n}A}(\vec{B})$$

$$\approx T_A(\vec{B}),$$

$$T_{AB}(\vec{B})$$

$$\approx {\int }_{-\infty }^{\infty }{\mathrm{d}}^2\vec{s}{T}_A(\vec{s}+\vec{B}/2)T_B(\vec{s}-\vec{B}/2),$$

$$T_A(\vec{s})\equiv {\int }_{-\infty }^{\infty }\mathrm{d}z{\rho }^A(\vec{s},z),$$

$$\int {\mathrm{d}}^2\vec{s}{T}_A(\vec{s})=A.$$

$${\sigma }_{AB}^{\mathrm{total}}=\int {\mathrm{d}}^2\vec{B}\sum _{k=1}^{AB}{\mathcal{P}}_k(\vec{B}),$$

$${\mathcal{P}}_k(\vec{B})$$

$$\times \sum _{{\alpha }_{11}=0}^1\dots \sum _{{\alpha }_{AB}=0}^1\left[p_{11}({\alpha }_{11})p_{12}({\alpha }_{12})\cdots p_{AB}({\alpha }_{AB})\right]{\delta }_{k,{\alpha }_{11}+\dots +{\alpha }_{AB}}.$$

$$p_{ij}({\alpha }_{ij})\equiv {\left(t_{ij}{\sigma }_{\mathrm{nn}}^{\mathrm{total}}\right)}^{{\alpha }_{ij}}{\left(1-t_{ij}{\sigma }_{\mathrm{nn}}^{\mathrm{total}}\right)}^{1-{\alpha }_{ij}}$$

$$t_{ij}\equiv t_{\mathrm{nn}}\left(\vec{B}+{\vec{S}}_i^A-{\vec{S}}_j^B\right).$$

$${\sigma }_{\mathrm{nn}}^{\mathrm{total}}=\int {\mathrm{d}}^2\vec{b}\sum _{k=1}^{\infty }{p}_k(\vec{b}),$$

$$p_k(\vec{b})=\exp \left[-t_{\mathrm{nn}}(\vec{b}){\sigma }_{\mathrm{nn}}\right]\frac{{\left[t_{\mathrm{nn}}(\vec{b}){\sigma }_{\mathrm{nn}}\right]}^k}{k!}.$$

$${\sigma }_{\mathrm{nn}}$$

$$\frac{\mathrm{d}{\sigma }_{AB\to a+b+X}}{\mathrm{d}{\vec{p}}^a\mathrm{d}{\vec{p}}^b}=\int \mathrm{d}\vec{B}\sum _{k=1}^{AB}{\mathcal{P}}_k(\vec{B})$$

$$\prod _{r=1}^k\int {\mathrm{d}}^2{\vec{b}}_r\sum _{k_r=1}^{\infty }\frac{p_{k_r}({\vec{b}}_r)}{{\sigma }_{\mathrm{nn}}^{\mathrm{total}}}\prod _{\ell =1}^{k_r}\left[\frac{1}{{\sigma }_{\mathrm{nn}}}\sum _{f_{r\ell }}\int \mathrm{d}{\mathrm{PS}}_{f_{r\ell }}\frac{\mathrm{d}{\sigma }_{nn\to f_{r\ell }}^{\mathrm{sps}}}{\mathrm{d}{\mathrm{PS}}_{f_{r\ell }}}\right]$$

$$\displaystyle\left[\sum_{i=1}^{k}\sum_{j=1}^{k_{r}}\sum_{n}\delta^{\,(3)}\left(\vec{p}^{\,a}-\vec{p}_{ij}^{\,a_{n}}\right)\times\left\{\begin{array}[]{c}1,\ {\rm if}\ a_{n}\in f_{ij}\\ 0,\ {\rm if}\ a_{n}\notin f_{ij}\end{array}\right\}\right]$$

$$\displaystyle\left[\sum_{i=1}^{k}\sum_{j=1}^{k_{r}}\sum_{n}\delta^{\,(3)}\left(\vec{p}^{\,b}-\vec{p}_{ij}^{\,b_{n}}\right)\times\left\{\begin{array}[]{c}1,\ {\rm if}\ b_{n}\in f_{ij}\\ 0,\ {\rm if}\ b_{n}\notin f_{ij}\end{array}\right\}\right]\,.$$

$${\mathcal{P}}_k(\vec{B})\times \left[\prod _{r=1}^k{p}_{k_r}({\vec{b}}_r)\right]\times \left[\prod _{\ell =1}^{k_r}\frac{1}{{\sigma }_{\mathrm{nn}}}\frac{\mathrm{d}{\sigma }_{nn\to f_{r\ell }}^{\mathrm{sps}}}{\mathrm{d}{\mathrm{PS}}_{f_{r\ell }}}\right],$$

$$\frac{\mathrm{d}{\sigma }_{nn\to a+X}^{\mathrm{sps}}}{{\mathrm{d}}^3{\vec{p}}^a}$$

$$\displaystyle\times\sum_{i}\delta^{\,(3)}\left(\vec{p}^{\,a}-\vec{p}^{\,a_{i}}\right)\times\left\{\begin{array}[]{c}1,\ {\rm if}\ a_{i}\in f\\ 0,\ {\rm if}\ a_{i}\notin f\end{array}\right\}\,,$$

$$\frac{\mathrm{d}{\sigma }_{nn\to a+b+X}^{\mathrm{sps}}}{{\mathrm{d}}^3{\vec{p}}^a{\mathrm{d}}^3{\vec{p}}^b}$$

$$\displaystyle\times\sum_{i}\delta^{\,(3)}\left(\vec{p}^{\,a}-\vec{p}^{\,a_{i}}\right)\times\left\{\begin{array}[]{c}1,\ {\rm if}\ a_{i}\in f\\ 0,\ {\rm if}\ a_{i}\notin f\end{array}\right\}$$

$$\displaystyle\times\sum_{i}\delta^{\,(3)}\left(\vec{p}^{\,b}-\vec{p}^{\,b_{i}}\right)\times\left\{\begin{array}[]{c}1,\ {\rm if}\ b_{i}\in f\\ 0,\ {\rm if}\ b_{i}\notin f\end{array}\right\}\,.$$

$$\mathrm{d}{\sigma }_{nn\to a+X}^{\mathrm{sps}}=\sum _{ijk}{f}_i({\mu }_{\mathrm{fact}}^2)$$

$$\otimes f_j({\mu }_{\mathrm{fact}}^2)\otimes D_{k\to a}({\mu }_{\mathrm{frag}}^2).$$

$$\mathrm{d}{\sigma }_{\mathrm{pp}\to a+b+X}^{\mathrm{sps}}$$

$$\otimes f_j({\mu }_{\mathrm{fact}}^2)\otimes D_{k\to a}({\mu }_{\mathrm{frag}}^2)\otimes D_{l\to b}({\mu }_{\mathrm{frag}}^2).$$

$$0.5\le \frac{{\mu }_{\mathrm{fact}}}{{\mu }_{\mathrm{ren}}}\le 2,0.5\le \frac{{\mu }_{\mathrm{frag}}}{{\mu }_{\mathrm{ren}}}\le 2,$$

$$D_{\mathrm{c}}^{{\mathrm{D}}^0}$$

Full text

††journal: Physics Letters B

1 Introduction

The recent measurements of inclusive open heavy-flavour – particularly D and B mesons – in proton-proton (p-p) collisions at the CERN Large Hadron Collider (LHC) [1, 2, 3, 4, 5, 6, 7, 8] provide opportunities to expose different facets of Quantum Chromodynamics (QCD) [9]. On one hand, due to the heavy-quark mass which serves as a hard interaction scale, the perturbative QCD calculations [10, 11, 12, 13, 14, 15] can be extended e.g. to very small transverse momenta ($p_{\rm T}$) where the calculations with massless quarks become inherently invalid. As the measurements at low $p_{\rm T}$ are statistically very precise, they offer an ideal testbed to benchmark perturbative calculations at low interaction scales. On the other hand, open heavy-flavour production can be used as a tool to probe non-perturbative aspects of heavy-quark fragmentation [16] and the quark-gluon structure of protons and nuclei [17, 18, 19, 20, 21, 22, 23, 24]. The low-$p_{\rm T}$ open heavy-flavour production in proton-lead (p-Pb) collisions [25, 26, 27, 28] may also open prospects to disentangle non-linear saturation [29, 30] and collinearly factorized QCD pictures in a regime where both should be valid descriptions.

The inclusive production of two D mesons provides exciting further opportunities. While the heavy-quarks are predominantly produced in pairs, the experimental overall reconstruction efficiencies for two D-meson final states are low and, roughly, only one out of million primarily produced double-D events can be reconstructed. Nevertheless, simultanous production of two D mesons has been observed in p-p [31] and p-$\overline{\rm p}$ [32] collisions. This offers possibilities to test e.g. the heavy-quark vs. heavy-antiquark asymmetries [33] and, in particular, to study the double-parton scattering (DPS) [34, 35, 36, 37, 38]. While the formal theory of factorization in DPS has recently advanced significantly [39, 40, 41], we still know relatively little of the non-perturbative structure of e.g. the double-parton distributions (dPDFs) [42] which would be required in precise phenomenological applications. Thus, simplifying assumptions concerning DPS have to be made which can lead to apparent shortcomings. For instance, measurements are often interpreted in terms of an effective cross section $\sigma_{\rm eff}$ whose inverse is proportional to the DPS probability. Its value has been observed to differ significantly depending from which observable it is extracted [43]. It is conceivable that this is due to overly simplifying the problem of DPS or, alternatively, overlooking the contributions of single-parton scattering (SPS) [44]. A complementary approach to hard DPS is provided by Monte-Carlo event generators in which the soft or semi-hard multiparton interactions (MPIs) give rise to the underlying event found necessary to describe the multiplicity distributions in hadronic high-energy collisions [41, 45].

The generic prediction is that in proton-nucleus (p-$A$) collisions the DPS signal gets enhanced in comparison to p-p case, due to the possibility of proton to interact with two or more nucleons simultaneously [46, 47, 48, 49, 50]. As in p-p collisions, multiple interactions are necessary to explain the multiplicity distributions in collisions involving heavy nuclei [51, 52], but a clean experimental confirmation for hard DPS processes is still lacking. As we will conclude later in this letter, it seems realistic that the double D-meson production could provide the first direct evidence of DPS in p-Pb collisions at clearly perturbative scales, and that the signal should be visible already in the collected Run-II data. In reaching this conclusion we have first confronted our QCD framework with the LHCb p-p data and a reasonble agreement there encourages us to apply it to p-Pb collisions. Before getting into the actual results we will first describe our theoretical framework in the next two sections.

2 Double-inclusive production in nuclear collisions

We will estimate the double-inclusive cross sections in collision of two nuclei, $A$ and $B$, in terms of inclusive per-nucleon SPS cross sections $\sigma^{\rm sps}_{\rm nn\rightarrow\mathcal{O}+X}$ as

[TABLE]

where $\vec{p}^{\,a}$ and $\vec{p}^{\,b}$ refer to the momenta of the produced particles $a$ and $b$. If $a$ and $b$ are identical particles $m=1/2$, and $m=1$ otherwise. In the case of independent partonic interactions, the effective cross section $\sigma_{\rm eff}^{AB}$ in $A$-$B$ collision is process independent and can be interpreted as a purely geometric object,

[TABLE]

Here,

[TABLE]

where $t_{\rm nn}(\vec{b})$ is the overlap function between two nucleons at fixed impact parameter $\vec{b}$. In geometric sense, we would write

[TABLE]

where $t_{\rm n}({\vec{s}})$ is the transverse profile of nucleons obtained by integrating the density of nucleons $\rho^{\rm n}$ over the longitudinal spatial component,

[TABLE]

The overlap functions $T_{\mathrm{n}A}(\vec{B})$ and $T_{AB}(\vec{B})$ at fixed impact parameter $\vec{B}$ are here defined as [53]

[TABLE]

where the approximation holds for point-like nucleons, and

[TABLE]

where $T_{A}(\vec{S})$ is the standard nuclear thickness function

[TABLE]

and $\rho^{A}$ denotes the density of nuclei. In our notation, the normalization is

[TABLE]

Typically, the DPS contribution in Eq. (1) is derived [47, 48, 49, 54] by writing the DPS cross section in terms of dPDFs, and assuming that the dPDFs factorize into a product of single-parton PDFs and that the partonic cross sections for the two subprocesses are unrelated and spatially independent. Alternatively, Eq. (1) can be derived from an eikonal model for multiparton interactions. Indeed, in a Glauber-type approach, the total cross section for a collision of nuclei A and B can be written in the impact-parameter space as

[TABLE]

where $\mathcal{P}_{k}(\vec{B})$ is the probability of exactly $k$ nucleon-nucleon interactions at fixed impact parameter $\vec{B}$,

[TABLE]

In this expression, we have defined

[TABLE]

where $t_{ij}$ is an abbreviation for the overlap function between two nucleons

[TABLE]

The second line in Eq. (11) thus corresponds to the probability of getting exactly $k$ nucleon-nucleon interactions (and $AB-k$ missing ones) at fixed geometric configuration. The total cross section in a single nucleon-nucleon collisions $\sigma_{\rm nn}^{\rm total}$ is given by

[TABLE]

where the probability $p_{k}(\vec{b})$ for $k$ partonic interactions is considered to be Poissonian,

[TABLE]

The quantity $\sigma_{\rm nn}$ appearing in Eq. (15) is the integrated inclusive cross sections,

[TABLE]

where the summation is over all exclusive final states $f$. We will always make a distinction between the (intensive) total cross section like $\sigma_{\rm nn}^{\rm total}$, and (extensive) integrated cross section like $\sigma_{\rm nn}$. The double-inclusive cross section can now be written as

[TABLE]

In the equation above, each term is a product of the form

[TABLE]

corresponding to the total probability density of having $k$ nucleon-nucleon interactions, each with exactly $k_{r=1,\ldots,k}$ partonic interactions resulting with a specific (exclusive) final state $f_{r\ell}$. The last two lines in Eq. (17) simply select those final states which contain the desired particles carrying the momenta $\vec{p}^{\,a}$ and $\vec{p}^{\,b}$, and the summation over $n$ accounts for the fact that the final state can contain several a or b particles. With some combinatorics, Eq. (17) simplifies to Eq. (1) when we identify

[TABLE]

and

[TABLE]

From the same formalism, also three-particle (in general $n$-particle) inclusive cross sections [55, 56, 49] can be derived.

3 Perturbative-QCD framework for open heavy flavour

In this paper we will be mostly concerned in the D-meson production at $p_{\rm T}>3\,{\rm GeV}$, which is the kinematic region considered in the LHCb double-D measurement [31]. In this region the inclusive production of D mesons can be reliably described within general-mass variable-flavour-number scheme (GM-VNFS). Schematically, the cross sections are convolutions of PDFs $f_{i}(x,\mu_{\rm fact}^{2})$, partonic cross sections $d\hat{\sigma}$, and fragmentation functions (FFs) $D_{k\rightarrow h}(z,\mu_{\rm frag}^{2})$,

[TABLE]

For single-inclusive D-meson production this has been considered at next-to-leading order (NLO) QCD first in Ref. [14] within the so-called SACOT scheme [57, 58]. In the SACOT scheme, the partonic cross sections for contributions in which the partonic subprocess is initiated by a charm quark or the fragmenting parton is a light one, are independent of the charm-quark mass $m_{\rm charm}$. This leads, in general, to diverging cross sections towards $p_{\rm T}\rightarrow 0$. In an alternative SACOT-$m_{\rm T}$ scheme [15] this unphysical behaviour is resolved by accounting for the underlying kinematic constraint of heavy-quark production. In this work we use the SACOT-$m_{\rm T}$ variant, albeit in the considered $p_{\rm T}>3\,{\rm GeV}$ region, both schemes should be equivalent within the scale uncertainties. Our default choice for factorization ($\mu_{\rm fact}$), fragmentation ($\mu_{\rm frag}$) and renormalization ($\mu_{\rm ren}$) scales is $\mu^{2}_{\rm fact}=\mu^{2}_{\rm frag}=\mu^{2}_{\rm ren}={p_{\rm T}^{2}+m_{\rm charm}^{2}}$, where $p_{\rm T}$ refers to the D-meson transverse momentum.

The SPS contribution in which the two D mesons, $h_{1}$ and $h_{2}$, are simultaneously produced is of the form,

[TABLE]

For this process, no GM-VFNS calculation is available. Thus, we will resort to the zero-mass approximation available in the NLO diphox [59] (v.1.2) code. Taking into account the large scale uncertainties, this approximation should be sufficiently precise in the considered $p_{\rm T}>3\,{\rm GeV}$ region. However, the kinematical cuts applied in the considered LHCb measurement [31] ($p_{\rm T}>3\,{\rm GeV}$ and $2<y<4$) include also a problematic configuration in which the two D mesons are collinear. In a full SACOT-$m_{\rm T}$ description this contribution would be finite, scaling as $\log(m_{\rm charm}^{2})$ where the remaining $\log(m_{\rm charm}^{2})$ terms would still need to be resummed via di-hadron FFs [60]. In a zero-mass calculation, however, the cross sections diverge in the collinear configuration. Here, as a proxy for the full SACOT-$m_{\rm T}$ treatment, we have regulated our calculations by imposing a physical cut $(\hat{p}_{1}+\hat{p}_{2})^{2}>4m_{\rm charm}^{2}$ for the fragmenting partons’ four momenta $\hat{p}_{1,2}$. Our central choice for the QCD scales here is the average $p_{\rm T}$ of the produced two D mesons.

The dominant uncertainty in our calculations comes from the unknown higher-order (NNLO and beyond) contributions. As usual, we estimate the potential size of these corrections by varying the QCD scales as

[TABLE]

around the central scale choices and finding the combinations that give the highest and lowest prediction for each considered observable. As default, we do the scale variations in sync for the two contributions in Eq. (1), the single-inclusive and double-inclusive SPS cross sections (17 scale configurations in total). We use NNPDF3.1pch PDFs [63] in which the intrinsic charm component is zero at the mass threshold $\mu_{\rm fact}=m_{\mathrm{c}}=1.51\,{\rm GeV}$. The fragmentation functions for D0 and D+ are taken from the KKKS08 analysis [61] (see Ref. [64] for a very recent alternative). The KKKS08 FFs have been fitted to $e^{+}e^{-}$ data from different experiments. We have checked that while the fits to BELLE [65] and OPAL [66] data give essentially equally good descriptions of the inclusive LHCb D0 and D+ cross sections at $\sqrt{s}=7\,{\rm TeV}$ [1], the FFs fitted to CLEO data [67] clearly overshoot the LHCb data at high $p_{\rm T}$. This is demonstrated in the upper panels of Figure 1. However, we have found that the D0-to-D+ ratios which are almost exclusively sensitive to the FFs are clearly best described by the OPAL variant, which also gives a better description than the BELLE FFs of the CMS midrapidity data [8] at very-high $p_{\rm T}$ [68]. Thus, in this paper, we adopt the OPAL FFs from the KKKS08 package. For D${}_{\mathrm{s}}^{\pm}$ and $\Lambda_{\mathrm{c}}^{\pm}$ FFs we use BKK05 [62] analysis. While the LHCb and ALICE single-inclusive D${}_{\rm s}^{\pm}$ data [1, 6] are well consistent with these FFs, the $\Lambda_{\mathrm{c}}^{\pm}$ data [1] are underestimated by the BKK05 $\Lambda_{\mathrm{c}}^{\pm}$ FFs. Our comparisons with the LHCb data on $D_{\mathrm{s}}^{\pm}$ and $\Lambda_{\mathrm{c}}^{\pm}$ are shown in the bottom panels of Figure 1.

The KKKS08 and BKK05 FFs do not discriminate between charge-conjugate states, but are given as a sum. (e.g. $D^{{\rm D}^{0}+\overline{\rm D^{0}}}_{i}$). In what follows, however, we will need the D-meson FFs one by one. Taking the D0 states here as an example, we will use the following prescription for the charm-quark containing state,

[TABLE]

and an analogous one for the antiquark-containing state,

[TABLE]

In addition to the NLO QCD framework described above, we present the predictions from Pythia 8 Monte-Carlo event generator using the standard “Monash 2013 tune” [69]. A sample of minimum-bias events were generated, including also MPIs, from which the different D-meson combinations within the LHCb acceptance were picked up to obtain the cross sections for each pair. In line with the LHCb measurements, each pair of D mesons is counted separately. In Figure 1 we also show the Pythia predictions for the inclusive D mesons and $\Lambda_{\mathrm{c}}^{\pm}$, generated with the provided Rivet analysis [70]. In general, the Pythia setup overpredicts the LHCb D-meson measurements, and the disagreement is stronger for D± and D${}^{\pm}_{\mathrm{s}}$ than for D0. A similar behaviour has been recently observed in the case of jets containing a ${\rm D}^{0}$ meson [71]. The measured $\Lambda_{\mathrm{c}}^{\pm}$ cross sections are, in turn, underestimated by Pythia. In the Monash tune the parameters related to charm fragmentation were constrained using a limited set of LEP data. Partly the interpretation of these data is hindered by the large feed-down from B-mesons. Furthermore, the data is not sensitive to $\mathrm{g}\rightarrow\mathrm{c\overline{c}}$ branchings that are abundant at the LHC. Thus the observed disagreement could potentially be cured by re-tuning the relevant parameters using a larger sample of charm-production data from LEP, HERA and LHC.

4 Results

We will now compare our results for double D-meson production with the LHCb p-p data [31], and make predictions for p-Pb collisions. As for $\sigma_{\rm eff}$, we will consider the variation $10~{}\text{mb}<\sigma_{\rm eff}<25~{}\text{mb}$ which is roughly the range deduced from jet, W± and photon measurements [43]. The uncertainty estimates shown in the plots combine the scale uncertainty and the variation in $\sigma_{\rm eff}$.

4.1 p-p collisions

In the case of p-p collision, Eqs. (1) and (2) reduce to

[TABLE]

Our results for the integrated cross sections within the LHCb acceptance are shown in Figure 2. For the opposite-sign D mesons (upper panel), the SPS contribution is clearly larger than the DPS one, and the agreement with the data is very good. The measured systematics among different combinations of D (and $\Lambda_{\rm c}$) species is well reproduced by the used set of FFs. As the cross section accumulates from the lower end of the considered $p_{\mathrm{T}}$ range, the scale uncertainty is sizable and dominates over the variation in $\sigma_{\rm eff}$. For like-sign final states (lower panel) the DPS becomes the dominant production mechanism. Again, the calculation agrees with the data within the scale uncertainties, though our central scale choice seem to somewhat overestimate the cross sections. The disagreement between the data and Pythia results is considerably larger than in the single-inclusive case (Figure 1). Apart from pairs including $\Lambda_{\rm c}$ the predicted cross sections are 4–8 times higher than the data. For pairs including $\Lambda_{\rm c}$ the systematic is again the opposite.

More insight can be obtained from Figure 3 where we show cross-section ratios. The upper panel shows ratios between the double like-sign vs. opposite-sign cross sections,

[TABLE]

These measure essentially the ratio between the DPS and SPS contributions. There is clearly a fair data-to-theory agreement within the scale and $\sigma_{\rm eff}$ uncertainties. Our central predictions somewhat overestimate the measured values which is consistent with Figure 2. The scale uncertainties do not cancel out since the partonic channels for like-sign and opposite-sign production are different (e.g. $\rm c\overline{\rm c}$ pair production is significant for $\rm D^{0}\overline{D^{0}}$ final state but not for $\rm D^{0}D^{0}$). Interestingly, the Pythia results are in excellent agreement with the LHCb data even though the absolute cross sections are way off. Since the numerator in Eq. (46) is sensitive to DPS (or MPIs in general), we conclude that the good agreement here suggest that the inconsistencies observed in Figure 2 are indeed due to poorly-constrained charm fragmentation, rather than the MPI modelling in Pythia [72, 73, 74].

The bottom panel of Figure 3 shows ratios

[TABLE]

From Eq. (45) we see that in the absence of SPS, this ratio would be equal to $\sigma_{\rm eff}$, but if there is a contribution from SPS, the ratio will be below $\sigma_{\rm eff}$. In general, our predictions for the opposite-sign case match very well with the data, but tend to underestimate the measured like-sign ratios. This is well in line with our earlier observations and also here a better overall agreement would be obtained if the DPS cross section would be somewhat smaller. Thus, the double-charm production data would prefer a somewhat larger phenomenological $\sigma_{\rm eff}$ than what other measurements indicate [43]. The Pythia predictions are here well compatible with our NLO calculations, though they somewhat undershoot the measured ratios both for like- and opposite-sign ratios. This further supports our conclusion that the disagreement observed in Figures 1 and 2 arise from the fragmentation scheme in Pythia.

4.2 p-Pb collisions

The reasonble description of the p-p data gives us confidence to apply the framework in p-$A$ collisions. In this case, Eqs. (1) and (2) reduce to

[TABLE]

with

[TABLE]

The impact-parameter integral for $A=208$ (Pb) gives

[TABLE]

taking $d=0.54\,{\rm fm}$ and $r=6.49\,{\rm fm}$ in the Woods-Saxon profile,

[TABLE]

and fixing $n_{0}$ by the normalization condition of Eq. (9). With $\sigma_{\rm eff}=10\ldots 25~{}{\rm mb}$, we find

[TABLE]

in full consistency e.g. with Ref. [47]. That is, the DPS signal is enhanced approximately by a factor of three in comparison to p-p scattering. Our results for the integrated cross sections within the LHCb kinematics are shown in Figure 4. Here, we have only considered D0 production which has the largest cross sections, see Figure 2, and the $y$ acceptance refers to that in the center-of-mass frame of the p-Pb collision. When computing the per-nucleon cross sections $\sigma^{\rm sps}_{nn\rightarrow a+b+X}$ and $\sigma^{\rm sps}_{nn\rightarrow a/b+X}$, we have used the EPPS16 nuclear modifications [75] for Pb. At the LHCb kinematics this leads to a $\sim 20\%$ suppression for p-Pb (forward) SPS cross sections, but since this is squared in DPS contribution, the suppression can reach $\sim 40\%$ in DPS case. For Pb-p configuration (backward) the nuclear-PDF effects are smaller. In comparison to the p-p case in Figure 2 the impact of enhanced DPS contribution is clear: Whereas in p-p case the DPS contribution to the opposite-sign yield was rather small in comparsion to SPS, in p-Pb collisions the two are comparable. For the like-sign yields the SPS contribution in p-Pb collisions is entirely overpowered by the DPS part, whereas in the p-p case the SPS still had a 20% contribution or so. Due to the additional contribution from the $T^{2}_{{\rm n}A}(\vec{B})$ integral in Eq. (49), the variation in $\sigma_{\rm eff}$ plays only a minor role as indicated in Figure 4. The $\sim$30% differences between forward and backward cross sections are due to the EPPS16 nuclear effects. Thus, by a suitable measurement where other theoretical uncertainties would cancel out, e.g. a forward-to-backward ratio for double D-meson production, further constraints for nuclear PDFs could, perhaps, be obtained.

An interesting question is whether these cross sections are large enough to be measured with the already collected Run-II data. In Run-II data taking the luminosities collected by the LHCb were $12.2\,{\rm nb}^{-1}$ for p-Pb (forward) and $18.6\,{\rm nb}^{-1}$ Pb-p (backward) collisions [27]. The overall detection efficiency $\epsilon$ for ${\rm D}^{0}\overline{{\rm D^{0}}}$ and ${\rm D}^{0}{{\rm D}}^{0}$ final states in the LHCb p-p measurement [31] was approximately $\epsilon\approx 1.2\times 10^{-6}$. Using these luminosities and efficiencies with our central theoretical predictions we calculate the expected number of events $N$ from which the statistical uncertainty is obtained as $\sqrt{N}/N$. These estimates are also shown in Figure 4. Within the scale uncertainties we expect approximately $10\dots 40$ ${\rm D}^{0}\overline{{\rm D}}^{0}$ pairs in p-Pb collisions (forward), and $20\ldots 80$ in Pb-p configuration (backward). For the like-sign case the corresponding numbers are $2\dots 20$ ${\rm D}^{0}{\rm D}^{0}$ pairs in p-Pb collisions (forward), and $4\ldots 40$ in Pb-p configuration (backward). Thus, we are led to conclude that the double D-meson production – at least the opposite-sign case – should be observable at the LHCb with the Run-II luminosity. Lowering the minimum-$p_{\rm T}$ cut below $3\,{\rm GeV}$ would easily increase the yields to a definitely measurable level, but towards lower $p_{\rm T}$ our predictions become increasingly uncertain.

As is well known, the increased importance of the DPS contribution in p-Pb collisions may significantly affect the kinematic distributions [76, 77]. Particularly interesting observable is the relative azimuthal-angle $\Delta\phi$ distribution of the two D mesons [78]. In p-p collisions [31] the $\Delta\phi$ distribution for ${\rm D}^{0}\overline{{\rm D^{0}}}$ peaks at $\Delta\phi=0$ for the logarithmically enhanced $g\rightarrow c\overline{c}$ splitting, and at $\Delta\phi=\pi$ due to the leading-order contributions that are back-to-back in transverse plane. These are commonly referred to as the near-side peak and the away-side peak, respectively. The disappearance of the away-side peak has long been predicted to be the smoking gun of saturation physics [79, 80]. However, the enhanced DPS contribution in p-Pb collisions will generate a $\Delta\phi$-independent contribution which levels off these peaks. Unfortunately our NLO QCD framework cannot reliably predict the $\Delta\phi$ dependence near $\Delta\phi=\pi$ but a soft-gluon resummation encoded e.g. in parton showers, would be required. To estimate the effect, we have fitted the $\Delta\phi$ dependence of the LHCb ${\rm D}^{0}\overline{{\rm D}^{0}}$ data [31] in p-p collisions assuming a negligible contribution from DPS. This assumption is consistent both with our results (see Figure 2) and also with Ref. [78] where it has been shown that for $0\lesssim\Delta\phi\lesssim\pi/2$ the fixed-order QCD quite correctly predicts the $\Delta\phi$ dependence. Our sketchy estimate for p-Pb is then a linear combination

[TABLE]

where the constant $\beta$ is determined by the relative importance of integrated SPS and DPS cross sections. Depending on the scale choices, we estimate the DPS contribution to be roughly between 20%…40%. The resulting projection for the $\Delta\phi$ dependence is shown in the lower panel of Figure 4 where the coloured band comes from the scale and $\sigma_{\rm eff}$ uncertainties. We see that both the near- and away-side peaks become less pronounced in p-Pb than what they are in p-p. Thus, we can confirm that the DPS contributions should be considered when interpreting the possible (probable?) weakening of the away-side two-particle correlations in terms of e.g. saturation physics. Since the DPS contribution should be nearly the same for ${\rm D}^{0}\overline{{\rm D}^{0}}$ and ${\rm D}^{0}{\rm D^{0}}$ final states (in our calculations they are equal), the difference

[TABLE]

should serve to subtract the “pedestal” DPS yield in a rather model-independent way and, as indicated above, correspond very closely to the SPS contribution in ${\rm D}^{0}\overline{{\rm D}^{0}}$ production. Even though the presented calculations are for double D-meson production, we would expect a similar reduction of the near- and away-side peaks also for light-flavour hadrons (such as $\pi^{+}\pi^{-}$) in due to enhanced DPS contribution in p-Pb collisions.

5 Summary

We have explored the double-inclusive D-meson production at the LHC with focus on the forward LHCb kinematics. The contributions of double-parton scatterings were included in the approximation of independent parton-parton collisions, and the required single-parton cross sections were computed within the collinearly factorized QCD at an NLO level. We confronted our predictions with the LHCb p-p data finding a good, or least an acceptable agreement within the QCD scale uncertainties and reasonable variation in the effective cross $\sigma_{\rm eff}$. As a whole, the LHCb data would prefer a rather large $\sigma_{\rm eff}$ compared to values derived from other final states. We also compared the LHCb p-p data with Pythia predictions. We found that the absolute cross sections for single- and double-inclusive open-charm production are not well reproduced by the widely used Monash tune. However, the cross-section ratios, which are less sensitive to the details in charm fragmentation, are described equally well or even better than what our NLO calculations do. Since the ratios are more sensitive to the multi-parton dynamics than the heavy-quark fragmentation model, it seems that the latter will need some further tuning to establish an agreement with the absolute cross sections.

In addition, we applied our framework to the case of p-Pb collisions in which the contribution from double-parton scattering is predicted to get significantly enhanced due to multiple nucleon-nucleon interactions. Our calculations, accounting for realistic reconstruction efficiencies, indicate that the yields should be high enough to be measured with the already-collected LHC Run-II data at the LHCb. This should provide a clear evidence for the hard double-parton scattering in p-Pb. As the contributions from single- and double-parton scatterings to opposite-sign double-D pair become comparable, also the azimuthal correlations are significantly altered. Therefore it seems necessary to take the double-parton scattering component into account when interpreting e.g. the possible complete or probable partial disappearance of the away-side peak.

Acknowledgments

We thank Michael Winn for clarifying us details of the LHCb p-p measurements and Peter Skands for discussions related to the Monash-tune applied in Pythia simulations. Our work was financed by the Academy of Finland, project n.o. 308301. The Finnish IT Center for Science (CSC) is acknowledged for computing resources within the project jyy2580 of T. Lappi.