Production of associated $\Upsilon$ and open charm hadrons in $pp$ collisions at $\sqrt{s}=7$ and $8$TeV via double parton scattering

Associated production of bottomonia and open charm hadrons in $pp$ collisions at $\sqrt{s}=7$ and $8$TeV is observed using data corresponding to an integrated luminosity of 3$fb^{-1}$ accumulated with the LHCb detector. The observation of five combinations, $\Upsilon(1S)D^0$, $\Upsilon(2S)D^0$, $\Upsilon(1S)D^+$, $\Upsilon(2S)D^+$ and $\Upsilon(1S)D^+_{s}$, is reported. Production cross-sections are measured for $\Upsilon(1S)D^0$ and $\Upsilon(1S)D^+$ pairs in the forward region. The measured cross-sections and the differential distributions indicate the dominance of double parton scattering as the main production mechanism. This allows a precise measurement of the effective cross-section for double parton scattering.


Introduction
Production of multiple heavy quark pairs in high-energy hadron collisions was first observed in 1982 by the NA3 collaboration in the channels π − (p) nucleon → J/ψ J/ψ + X [1,2]. Soon after, evidence for the associated production of four open charm particles in pion-nucleon reactions was obtained by the WA75 collaboration [3]. A measurement of J/ψ pair production in proton-proton (pp) collisions at √ s = 7 TeV [4] has been made by the LHCb collaboration in 2011. This measurement appears to be in good agreement with two models within the single parton scattering (SPS) mechanism, namely non-relativistic quantum chromodynamics (NRQCD) calculations [5] and k T -factorization [6]. However the obtained result also agrees with predictions [7] of the double parton scattering (DPS) mechanism [8][9][10][11][12].
The production of J/ψ pairs has also been observed by the D0 [13] and CMS [14] collaborations. A large double charm production cross-section involving open charm in pp collisions at √ s = 7 TeV has been observed by the LHCb collaboration [15]. The measured cross-sections exceed the SPS expectations significantly [16][17][18][19][20] and agree with the DPS estimates. A study of differential distributions supports a large role for the DPS mechanism in multiple production of heavy quarks.
The study of (bb)(cc) production in hadronic collisions started with the observation of B + c mesons in pp collisions by the CDF collaboration [21]. A detailed study of B + c production spectra in pp collisions by the LHCb collaboration [22] showed good agreement with leadingorder NRQCD calculations [23][24][25] including the SPS contribution only.
The leading-order NRQCD calculations using the same matrix element as in Ref. [23], applied to another class of (bb)(cc) production, namely associated production of bottomonia and open charm hadrons in the forward region, defined in terms of the rapidity y as 2 < y < 4.5, predict [26] R SPS = σ Υcc σ Υ = (0.2 − 0.6) % , where σ Υcc denotes the production cross-section for associated production of Υcc-pair and σ Υ denotes the inclusive production cross-section of Υ mesons. A slightly smaller value of R SPS is obtained through the k T -factorization approach [17,[27][28][29][30][31] using the transverse momentum dependent gluon density from Refs. [32][33][34], Within the DPS mechanism, the Υ meson and cc-pair are produced independently in different partonic interactions. Neglecting the parton correlations in the proton, the contribution of this mechanism is estimated according to the formula [35][36][37] where σ cc and σ Υ are the inclusive charm and Υ cross-sections, and σ eff is an effective cross-section, which provides the proper normalization of the DPS cross-section estimate.
Here we report the first observation of associated production of bottomonia and open charm hadrons. The production cross-sections and the differential distributions are measured. The latter provide crucial information for understanding the production mechanism. The analysis is performed using the Run 1 data set recorded by the LHCb detector, consisting of 1 fb −1 of integrated luminosity accumulated at √ s = 7 TeV and 2 fb −1 accumulated at 8 TeV.

Detector and data sample
The LHCb detector [43,44] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of siliconstrip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex, the impact parameter, is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger [45], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware stage, events for this analysis are selected requiring dimuon candidates with a product of their transverse momenta p T larger than 1.7 (2.6) GeV 2 /c 2 for data collected at √ s = 7 (8) TeV. In the subsequent software trigger, two well reconstructed tracks are required to have hits in the muon system, to have p T > 500 MeV/c and p > 6 GeV/c and to form a common vertex. Only events with a dimuon candidate with a mass m µ + µ − larger than 4.7 GeV/c 2 are retained for further analysis.
The simulation is performed using the LHCb configuration [46] of the Pythia 6 event generator [47]. Decays of hadronic particles are described by EvtGen [48] in which final-state photons are generated using Photos [49]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [50,51] as described in Ref. [52].

Event selection
The event selection strategy is based on the independent selection of Υ(1S), Υ(2S) and Υ(3S) mesons (jointly referred to by the symbol Υ throughout the paper) and charmed hadrons, namely D 0 , D + and D + s mesons and Λ + c baryons (jointly referred to by the symbol C herafter) originating from the same pp collision vertex. The Υ candidates are reconstructed via their dimuon decays, and the D 0 → K − π + , D + → K − π + π + , D + s → K + K − π + and Λ + c → pK − π + decay modes are used for the reconstruction of charm hadrons. Charge conjugate processes are implied throughout the paper. The fiducial region for this analysis is defined in terms of the p T and the rapidity y of Υ and C hadrons to be p Υ T < 15 GeV/c, 2.0 < y Υ < 4.5, 1 < p C T < 20 GeV/c and 2.0 < y C < 4.5. The event selection for Υ → µ + µ − candidates follows previous LHCb studies [41], and the selection of C hadrons follows Refs. [15,53]. Only good quality tracks [54], identified as muons [55], kaons, pions or protons [56] are used in the analysis. A good quality vertex is required for Υ → µ + µ − , D 0 → K − π + , D + → K − π + π + , D + s → K + K − π + and Λ + c → pK − π + candidates. For D + s → K + K − π + candidates, the mass of the K + K − pair is required to be in the region m K + K − < 1.04 GeV/c 2 , which is dominated by the D + s → φπ + decay. To suppress combinatorial background the decay time of C hadrons is required to exceed 100 µm/c. Full decay chain fits are applied separately for selected Υ and C candidates [57]. For Υ mesons it is required that the vertex is compatible with one of the reconstructed pp collision vertices. In the case of long-lived charm hadrons, the momentum direction is required to be consistent with the flight direction calculated from the locations of the primary and secondary vertices. The reduced χ 2 of these fits, both χ 2 fit (Υ) /ndf and χ 2 fit (C) /ndf, are required to be less than 5, where ndf is the number of degrees of freedom in the fit. The requirements favour the selection of charm hadrons produced promptly at the pp collision vertex and significantly suppress the feed down from charm hadrons produced in decays of beauty hadrons. The contamination of such C hadrons in the selected sample varies between (0.4 ± 0.2)% for D + mesons to (1.5 ± 0.5)% for Λ + c baryons. The selected Υ and C candidates are paired to form ΥC candidates. A global fit to the ΥC candidates is performed [57], similar to that described above, which requires both hadrons to be consistent with originating from a common vertex. The reduced χ 2 of this fit, χ 2 fit (ΥC) /ndf, is required to be less than 5. This reduces the background from the pile-up of two independent pp interactions producing separately a Υ meson and C hadron to a negligible level, keeping 100% of the signal Υ mesons and C hadrons from the same primary vertex. The two-dimensional mass distributions for ΥC pairs after the selection are displayed in Fig. 1.

Signal extraction and cross-section determination
The event yields are determined using unbinned extended maximum likelihood fits to the two-dimensional ΥC mass distributions of the selected candidates. The fit model is a sum of several components, each of which is the product of a dimuon mass distribution, corresponding to an individual Υ state or combinatorial background, and a C candidate mass distribution, corresponding to a C signal or combinatorial background component. The Υ(1S) → µ + µ − , Υ(2S) → µ + µ − and Υ(3S) → µ + µ − signals are each modelled by a double-sided Crystal Ball function [4,58,59] and referred to as S Υ in this section. A modified Novosibirsk function [60] (referred to as S C ) is used to describe the D 0 → K − π + , D + → K − π + π + , D + s → K + K − π + and Λ + c → pK − π + signals. All shape parameters and signal peak positions are fixed from fits to large inclusive Υ → µ + µ − and C hadron data samples. Combinatorial background components B µ + µ − and B C are modelled with a product of exponential and polynomial functions with a slope parameter β and a polynomial function P n , which is represented as a Bézier sum of basic Bernstein polynomials of order n with non-negative coefficients [61]. For the large yield ΥD 0 and ΥD + samples, the second-order polynomials (n = 2) are used in the fit, while n = 1 is used for the ΥD + s and ΥΛ + c cases.  These basic functions are used to build the components of the two dimensional mass fit following Ref. [15]. For each C hadron the reconstructed signal sample consists of the following components: -Three ΥC signal components: each is modelled by a product of the individual signal Υ components, -Three components describing the production of single Υ mesons together with combinatorial background for the C signal: each component is modelled by a product of the signal Υ component, S Υ (m µ + µ − ) and the background component B C (m C ).
-Single production of C hadrons together with combinatorial background for the Υ component: this is modelled by a product of the signal C component, S C (m C ), and the background component B µ + µ − (m µ + µ − ). -Combinatorial background: this is modelled by a product of the individual background components B µ + µ − (m µ + µ − ) and B C (m C ).
For each C hadron the complete fit function F(m µ + µ − , m C ) is where the different coefficients N ΥC , N ΥB , N BC and N BB are the yields of the eight components described above.
The fit results are summarized in Table 1, and the fit projections are presented in Figs. 2, 3, 4 and 5. The statistical significances of the signal components are determined using a Monte-Carlo technique with a large number of pseudoexperiments. They are presented in Table 2. For the five modes, Υ(1S)D 0 , Υ(2S)D 0 , Υ(1S)D + , Υ(2S)D + and Υ(1S)D + s , the significances exceed five standard deviations. No significant signals are found for the associated production of Υ mesons and Λ + c baryons. The possible contribution from pile-up events is estimated from data following the method from Refs. [15,53] by relaxing the requirement on χ 2 fit (ΥC) /ndf. Due to the requirements χ 2 fit (Υ) /ndf < 5 and χ 2 fit (C) /ndf < 5, the value of χ 2 fit (ΥC) /ndf does not exceed 5 units for signal events with Υ and C hadron from the same pp collision vertex. The background is subtracted using the sPlot technique [63]. The χ 2 fit (ΥC) /ndf distributions are shown in Fig. 6. The distributions exhibit two components: the peak at low χ 2 is attributed to associated ΥC production, and the broad structure at large values of χ 2 corresponds to the contribution from pile-up events. The distributions are fitted with a function that has two components, each described by a Γ-distribution. The shape is Table 2: Statistical significances of the observed ΥC signals in units of standard deviations determined using pseudoexperiments. The values in parentheses indicate the statistical significance calculated using Wilks' theorem [62].
motivated by the observation that χ 2 fit /ndf should follow a scaled-χ 2 distribution. The possible contribution from pile-up events is estimated by integrating the pile-up component in the region χ 2 fit (ΥC) /ndf < 5. It does not exceed 1.5% for all four cases and is neglected. The production cross-section is determined for the four modes with the largest yield: Υ(1S)D 0 , Υ(2S)D 0 , Υ(1S)D + and Υ(2S)D + . The cross-section is calculated using a subsample of events where the reconstructed Υ candidate is explicitly matched to the dimuon candidate that triggers the event. This requirement reduces signal yields by approximately 20%, but allows a robust determination of trigger efficiencies. The cross-section for the associated production of Υ mesons with C hadrons in the kinematic range of LHCb is calculated as where L is the integrated luminosity [64], B µ + µ − and B C are the world average branching fractions of Υ → µ + µ − and the charm decay modes [42], and N ΥC corr is the efficiencycorrected yield of the signal ΥC events in the kinematic range of this analysis. Production cross-sections are determined separately for data sets accumulated at √ s = 7 and 8 TeV. The efficiency-corrected signal yields N ΥC corr are determined using an extended unbinned maximum likelihood fit to the weighted two-dimensional invariant mass distributions of the selected ΥC candidates. The weight ω for each event is calculated as ω = 1/ε tot , where ε tot is the total efficiency for the given event.
The effective DPS cross-section and the ratios R ΥC are calculated as where σ Υ is the production cross-section of Υ mesons taken from Ref. [41]. The double-differential production cross-sections of charm mesons has been measured at √ s = 7 TeV in the region 2.0 < y C < 4.5, p C T < 8 GeV/c [38]. According to FONLL calculations [65][66][67], the contribution from the region 8 < p C T < 20 GeV/c is significantly smaller than the uncertainty for the measured cross-section in the region 1 < p C T < 8 GeV/c. It allows to estimate the production cross-section of charm mesons in the region 2.0 < y C < 4.5, Projections from two-dimensional extended unbinned maximum likelihood fits in bands a) 1.844 < m K − π + < 1.887 MeV/c 2 , b) 9.332 < m µ + µ − < 9.575 GeV/c 2 , c) 9.889 < m µ + µ − < 10.145 GeV/c 2 and d) 10.216 < m µ + µ − < 10.481 GeV/c 2 . The total fit function is shown by a solid thick (red) curve; three individual ΥD 0 signal components are shown by solid thin (red) curves; three components describing Υ signals and combinatorial background in K − π + mass are shown with short-dashed (blue) curves; the component modelling the true D 0 signal and combinatorial background in µ + µ − mass is shown with a long-dashed (green) curve and the component describing combinatorial background is shown with a thin dotted (black) line.
. For the production cross-section of charm mesons at √ s = 8 TeV, the measured cross-section at √ s = 7 TeV is rescaled by the ratio R FONLL 8/7 (p T , y) of the double-differential cross-sections, as calculated with FONLL [65][66][67] at √ s = 8 and 7 TeV.  Figure 3: Projections from two-dimensional extended unbinned maximum likelihood fits in bands a) 1.848 < m K − π + π + < 1.891 MeV/c 2 , b) 9.332 < m µ + µ − < 9.575 GeV/c 2 , c) 9.889 < m µ + µ − < 10.145 GeV/c 2 and d) 10.216 < m µ + µ − < 10.481 GeV/c 2 . The total fit function is shown by a solid thick (red) curve; three individual ΥD + signal components are shown by solid thin (red) curves; three components describing Υ signals and combinatorial background in K − π + π + mass are shown with short-dashed (blue) curves; the component modelling the true D + signal and combinatorial background in µ + µ − mass is shown with a long-dashed (green) curve and the component describing combinatorial background is shown with a thin dotted (black) line.
The ratios , defined in Eq. (6), are calculated as where ε ΥC denotes the average efficiency. Within the DPS mechanism, the transverse momenta and rapidity spectra of C mesons for the signal Υ(1S)C and Υ(2S)C events are expected to be the same. This allows to express the ratio of the average ε ΥC efficiencies  Figure 4: Projections from two-dimensional extended unbinned maximum likelihood fits in bands a) 1.952 < m (K − K + ) φ π + < 1.988 MeV/c 2 , b) 9.332 < m µ + µ − < 9.575 GeV/c 2 , c) 9.889 < m µ + µ − < 10.145 GeV/c 2 and d) 10.216 < m µ + µ − < 10.481 GeV/c 2 . The total fit function is shown by a solid thick (red) curve; three individual ΥD + s signal components are shown by solid thin (red) curves; three components describing Υ signals and combinatorial background in (K − K + ) φ π + mass are shown with short-dashed (blue) curves; the component modelling the true D + s signal and combinatorial background in µ + µ − mass is shown with a long-dashed (green) curve and the component describing combinatorial background is shown with a thin dotted (black) line.
in terms of ratio of average efficiencies for inclusive Υ mesons and the latter is taken from Ref. [41]. The total efficiency ε tot , for each ΥC candidate is calculated following Ref. [15] as and applied individually an on event-by-event basis, where ε tot Υ and ε tot C are the total  Figure 5: Projections from two-dimensional extended unbinned maximum likelihood fits in bands a) 2.273 < m pK − π + < 2.304 MeV/c 2 , b) 9.332 < m µ + µ − < 9.575 GeV/c 2 , c) 9.889 < m µ + µ − < 10.145 GeV/c 2 and d) 10.216 < m µ + µ − < 10.481 GeV/c 2 . The total fit function is shown by a solid thick (red) curve; three individual ΥΛ + c signal components are shown by solid thin (red) curves; three components describing Υ signals and combinatorial background in pK − π + mass are shown with short-dashed (blue) curves; the component modelling the true Λ + c signal and combinatorial background in µ + µ − mass is shown with a long-dashed green curve and the component describing combinatorial background is shown with a thin dotted (black) line. efficiencies for Υ and charm hadrons respectively. These efficiencies are calculated as where ε rec is the detector acceptance, reconstruction and event selection efficiency and ε trg is the trigger efficiency for selected events. The particle identification efficiencies for where ε ID µ ± , ε ID K and ε ID π are the efficiencies for the single muon, kaon and pion identification, respectively.
The efficiencies ε rec and ε trg are determined using simulated samples of Υ, D 0 and D + events as a function of p T and y of the Υ and the C hadron. The differential treatment results in a robust determination of the efficiency-corrected signal yields, with no dependence on the particle spectra in the simulated samples. The derived values of the efficiencies are corrected to account for small discrepancies in the detector response between data and simulation. These corrections are obtained using data-driven techniques [54,55].
The efficiencies for muon, kaon and pion identification are determined directly from data using large samples of low-background J/ψ → µ + µ − and D * + → (D 0 → K − π + ) π + decays. The identification efficiencies are evaluated as a function of the kinematic parameters of the final-state particles, and the track multiplicity in the event [56].
The efficiency is dependent on the polarisation of the Υ mesons [41,59,68,69] The polarisation of the Υ mesons produced in pp collisions at √ s = 7 TeV at high p Υ T and central rapidity has been studied by the CMS collaboration [70] in the centre-of-mass helicity, Collins-Soper [71] and the perpendicular helicity frames. No evidence of significant transverse or longitudinal polarisation has been observed for the region 10 < p Υ T < 50 GeV/c, y Υ < 1.2. Therefore, the efficiencies are calculated under the assumption of unpolarised production of Υ mesons and no corresponding systematic uncertainty is assigned on the cross-section. Under the assumption of transversely polarised Υ mesons with λ ϑ = 0.2 in the LHCb kinematic region, 1 the total efficiency would result in an decrease of 3% [41].

Kinematic distributions of ΥC events
The differential distributions are important for the determination of the production mechanism. In this section, the shapes of differential distributions for Υ(1S)D 0 and Υ(1S)D + events are studied. Assuming that the production mechanism of ΥC events is essentially the same at √ s = 7 and 8 TeV, both samples are treated together in this section. The normalized differential distribution for each variable v is calculated as where N ΥC corr,i is the number of efficiency-corrected signal events in bin i of width ∆v, and N ΥC corr is the total number of efficiency-corrected events. The differential distributions are presented for the following variables  . The transverse momentum spectra, derived within the DPS mechanism using the measurements from Refs. [38,41], are shown with the open (blue) squares. The SPS predictions [72] for the p Υ T spectra are shown with dashed (orange) and long-dashed (magenta) curves for calculations based on the k T -factorization and the collinear approximation, respectively. All distributions are normalized to unity.
The distributions are shown in Figs. 7, 8, 9, 10 and 11. Only statistical uncertainties are displayed on these figures, as the systematic uncertainities discussed in Sect. 6 are small. For all variables the width of the resolution function is much smaller than the bin width, i.e. the results are not affected by bin-to-bin migration.  Figure 8: Background-subtracted and efficiency-corrected y Υ (top) and y C (bottom) distributions for Υ(1S)D 0 (left) and Υ(1S)D + (right) events. The rapidity spectra, derived within the DPS mechanism using the measurements from Refs. [38,41], are shown with the open (blue) squares. The SPS predictions [72] for the y Υ spectra are shown with dashed (orange) and long-dashed (magenta) curves for calculations based on the k T -factorization and the collinear approximation, respectively. All distributions are normalized to unity.
The shapes of the measured differential distributions are compared with the SPS and DPS predictions. The DPS predictions are deduced from the measurements given in Refs. [38,41], using a simplified simulation assuming uncorrelated production of the Υ(1S) and charm hadron. The agreement between all measured distributions and the DPS predictions is good. For the SPS mechanism, the predictions [72] based on k T -factorization [17,[27][28][29][30][31] using the transverse momentum dependent gluon density from Refs. [32][33][34] are used along with the collinear approximation [26] with the leading-order gluon density taken from Ref. [73]. The transverse momentum and rapidity distributions of Υ(1S) mesons also agree well with SPS predictions based on k T -factorization, while the shape of the transverse momentum spectra of Υ mesons disfavours the SPS predictions obtained using the collinear approximation. The shapes of the y Υ distribution have very  [72] for the shapes of ∆φ distribution are shown with dashed (orange) and long-dashed (magenta) curves for calculations based on the k T -factorization and the collinear approximation, respectively. The solid (blue) curves in the ∆y plots show the spectra obtained using a simplified simulation based on data from Refs. [38,41]. The dashed (green) lines show the triangle function expected for totally uncorrelated production of two particles, uniformly distributed in rapidity. All distributions are normalized to unity.
limited sensitivity to the underlying production mechanism. The distribution |∆φ| is presented in Fig. 9(a,b). The DPS mechanism predicts a flat distribution in ∆φ, while for SPS a prominent enhancement at |∆φ| ∼ π is expected in collinear approximation. The enhancement is partly reduced taking into account transverse momenta of collinding partons [30,74] and it is expected to be further smeared out at next-to-leading order. The measured distributions for Υ(1S)D 0 and Υ(1S)D + events, shown in Fig. 9(a,b) agree with a flat distribution. The fit result with a constant function gives a p-value of 6% (12%) for the Υ(1S)D 0 (Υ(1S)D + ) case, indicating that the SPS contribution to the data is small. The shape of ∆y distribution is defined primarily by the acceptance of LHCb experiment 2 < y < 4.5 and has no sensitivity to the underlying production mechanism, in the limit of current statistics.

Systematic uncertainties
The systematic uncertainties related to the measurement of the production cross-section for ΥC pairs are summarized in Table 3 and discussed in detail in the following.
The signal shapes and parameters are taken from fits to large low-background inclusive Υ → µ + µ − and charm samples. The parameters, signal peak positions and resolutions and the tail parameters for the double-sided Crystal Ball and the modified Novosibirsk functions, are varied within their uncertainties as determined from the calibration samples. The small difference in parameters between the data sets obtained at √ s = 7 and 8 TeV is also used to assign the systematic uncertainty. For D 0 and D + signal peaks alternative fit models have been used, namely a double-sided asymmetric variant of an Apolonious function [75] without power-law tail, a double-sided Crystal Ball function and an asymmetric Student-t shape. The systematic uncertainty related to the parameterization of the combinatorial background is determined by varying the order of the polynomial function in Eq. (7) between zeroth and second order. For the purely combinatorial background component (last line in Eq. (8)), a non-factorizable function is used where the parameters β 1 , β 2 and κ i,j are allowed to float in the fit, and P i n and P j k are basic Bernstein polynomials, and the order of these polynomials, n and k, is varied between zero Table 3: Summary of relative systematic uncertainties for σ ΥC (in %). The total systematic uncertainty does not include the systematic uncertainty related to the knowledge of integrated luminosity [64]. The symbol ⊕ denotes the sum in quadrature.
Source and two. The corresponding variations of ΥC signal yields are taken as the systematic uncertainty related to the description of the signal and background components.
Other systematic uncertainties are related to the imperfection of the Photos generator [49] to describe the radiative tails in Υ → µ + µ − decays. This systematic is studied in Ref. [76] and taken to be 1%.
The systematic uncertainty related to efficiency correction is estimated using an alternative technique for the determination of N ΥC corr , where the efficiency-corrected yields are obtained via where w i is the signal event weight, obtained with the sPlot technique [63] using fits to the efficiency-uncorrected data sets, and ε tot is a total efficiency for the given event, defined with Eq. (13). The difference in the efficiency-corrected yields with respect to the baseline approach of 0.1 (1.3)% for ΥD 0 (ΥD + ), is assigned as the corresponding systematic uncertainty. The systematic uncertainty related to the particle identification is estimated to be 0.2% for muons and 0.5 (0.8)% for hadrons for the Υ(1S)D 0 (Υ(1S)D + ) case and is obtained from the uncertainties for the single particle identification efficiencies using an error propagation technique with a large number of pseudoexperiments. The same approach is Table 4: Summary of relative systematic uncertainties for σ eff (in %). The reduced uncertainty for C hadron production cross-section, denoted as δ(σ C ), is recalculated from Ref. [38] taking into account the cancellation of correlated systematic uncertainties.  Table 5: Summary of relative systematic uncertainties for the ratios R ΥC and R D 0 /D + (in %).
Source used to propagate the uncertainties in ε acc , ε rec and ε trg related to the limited simulation sample size. The efficiency is corrected using data-driven techniques to account for small differences in the tracking efficiency between data and simulation [54,55]. The uncertainty in the correction factor is propagated to the cross-section measurement using pseudoexperiments resulting in a global 0.4 (0.5)% systematic uncertainty for the ΥD 0 (ΥD + ) cases plus an additional uncertainty of 0.4% per track. The knowledge of the hadronic interaction length of the detector results in an uncertainty of 1.4% per final-state hadron [54].
The systematic uncertainty associated with the trigger requirements is assessed by studying the performance of the dimuon trigger for Υ(1S) events selected using the single muon high-p T trigger [45] in data and simulation. The comparison is performed in bins of the Υ(1S) meson transverse momentum and rapidity and the largest observed difference of 2.0% is assigned as the systematic uncertainty associated with the imperfection of trigger simulation [41].
Using large samples of low-background inclusive Υ → µ + µ − , D 0 → K − π + and D + → K − π + π + events, good agreement between data and simulation is observed for the selection variables used in this analysis, in particular for dimuon and charm vertex quality and χ 2 fit (Υ)/ndf. The small differences seen would affect the efficiencies by less than 1.0%, which is conservatively taken as a systematic uncertainty accounting for the disagreement between data and simulation.
The systematic uncertainty related to the uncertainties of the branching fractions of D 0 and D + mesons is 1.3% and 2.1% [42]. The integrated luminosity is measured using a beam-gas imaging method [77,78]. The absolute luminosity scale is determined with 1.7 (1.2)% uncertainty for the sample collected at √ s = 7 (8) TeV, dominated by the vertex resolution for beam-gas interactions, the spread of the measurements and the detector alignment [64,78,79].
The selection criteria favour the selection of charm hadrons produced promptly at the pp collision vertex and significantly suppress the feed down from charm hadrons produced in decays of beauty hadrons. The remaining feed down is estimated separately for DPS and SPS processes with the simultaneous production of an Υ meson and a bbpair. The former is estimated using simulation, normalized to the measured bb and cc production cross-sections [38,80] and validated using a data-driven technique. It is found to be smaller than 1.5% of the observed signal and is neglected. The contribution from SPS processes with the associated production of Υ meson and bb pairs is estimated using the prediction for the ratio of production cross-sections, obtained using the k T -factorization approach with the transverse momentum dependent gluon density taken from Refs. [32][33][34]. The uncertainty reflects the variation of scale and the difference with results obtained using the collinear approximation with the gluon density from Ref. [73]. Combining the estimates from Eqs. (19), (1) and (2) with the probability for a charm hadron from the decay of beauty hadron to pass the selection criteria, this feed down is found to be totally negligible. The effect of possible extreme polarization scenarios for Υ mesons from SPS processes is proportional to the SPS contamination, α SPS , and could lead to +0.08 (−0.16)α SPS correction [68] to the cross-sections σ ΥC and the ratios R ΥC for totally transverse (longitudinal) polarizations of Υ mesons in centre-of-mass helicity frame. It is very small for small SPS contamination. The corresponding corrections to ratios R D 0 /D + are non-zero only if SPS has different contributions to ΥD 0 and ΥD + production processes and accouts for +0.08 (−0.16)∆α SPS , where ∆α SPS is the difference in SPS contaminations to the considered processes. The same estimate is valid also for the ratios R Υ(1S)/Υ(2S) .
A large part of the systematic uncertainties cancels in the ratio R ΥC and in the variable σ eff . The systematic uncertainties for σ eff , R ΥC and R D 0 /D + are summarized in Tables 4 and 5. For the production cross-section of charm mesons at √ s = 8 TeV the measured cross-section at √ s = 7 TeV is extrapolated using FONLL calculations [65][66][67]. The uncertainty related to the imperfection of the extrapolation is estimated from the comparison of the measured ratio σ C √ s=13 TeV /σ C √ s=7 TeV [38,81] and the corresponding FONLL estimate. As a result of this comparison the C hadron production cross-section is scaled up by 2.7% and a systematic uncertainty of 2.1% is assigned. The systematic uncertainty for the ratios R is small compared to the statistical uncertainty and is neglected.

Results and discussion
The associated production of Υ and charm mesons is studied. Pair production of Υ(1S)D 0 , Υ(2S)D 0 , Υ(1S)D + , Υ(2S)D + and Υ(1S)D + s states is observed with significances exceeding five standard deviations. The production cross-sections in the fiducial region 2.0 < y Υ < 4.5, p Υ T < 15 GeV/c, 2.0 < y C < 4.5 and 1 < p C T < 20 GeV/c are measured for Υ(1S)D 0 and Υ(1S)D + final states at √ s = 7 and 8 TeV as: where the first uncertainty is statistical, and the second is the systematic uncertainty from Table 3, combined with the uncertainty related to the knowledge of the luminosity. All these measurements are statistically limited. The measured cross-sections are in agreement with the DPS expectations from Eq. (5), and significantly exceed the expectations from the SPS mechanism in Eqs. (1) and (2). Differential kinematic distributions are studied for ΥD 0 and ΥD + final states. All of them are in good agreement with DPS expectations as the main production mechanism.
The ratios of the cross-sections for Υ(1S)D 0 and Υ(1S)D + are where the systematic uncertainty is discussed in detail in Sect. 6. The results are compatible with the DPS expectation of 2.41 ± 0.18 from Eq. (6a). The cross-section ratios R ΥC are measured to be Extrapolating the ratios R ΥC down to p C T = 0 using the measured transverse momentum spectra of D 0 and D + mesons from Ref. [38], and using the fragmentation fractions f (c → D 0 ) = 0.565 ± 0.032 and f (c → D + ) = 0.246 ± 0.020, measured at e + e − colliders operating at a centre-of-mass energy close to the Υ(4S) resonance [82], the ratios R Υcc are calculated to be which significantly exceed SPS expectations from Eqs. (1) and (2).
The large statistical uncertainty for the other ΥC modes does not allow to obtain a numerical model-independent measurement, but, assuming similar kinematics for Υ(2S) and charm mesons to the prompt production, the following ratios are measured where B 2/1 is the ratio of dimuon branching fractions of Υ(2S) and Υ(1S) mesons and where the systematic uncertainties are negligible compared to statistical uncertainties. These values are smaller than, but compatible with the DPS expectations from Eq. (6b). For the ΥD + production one obtains where again the systematic uncertainties are negligible with respect to the statistical ones and are ignored. These values are compatible with the DPS expectation of 25% from Eq. (6b). Neglecting the contributions from SPS mechanism, the effective cross-section σ eff is determined using Eq. (10a) for the √ s = 7 TeV data as The mean value of is in good agreement with those obtained for √ s = 7 TeV data. Averaging these values, σ eff is found to be The large value of the cross-section for the associated production of Υ and open charm hadrons, compatible with the DPS estimate of Eq. (4), has important consequences for the interpretation of heavy-flavor production measurements, especially inclusive measurements and possibly for b-flavor tagging [89][90][91][92], where the production of uncorrelated charm hadrons could affect the right assignment of the initial flavour of the studied beauty hadron.

Summary
The associated production of Υ mesons with open charm hadrons is observed in pp collisions at centre-of-mass energies of 7 and 8 TeV using data samples corresponding to integrated luminosities of 1 fb −1 and 2 fb −1 respectively, collected with the LHCb detector. The production of Υ(1S)D 0 , Υ(2S)D 0 , Υ(1S)D + , Υ(2S)D + and Υ(1S)D + s pairs is observed with significances larger than 5 standard deviations. The production cross-sections in the fiducial region 2.0 < y Υ < 4.5, p Υ T < 15 GeV/c, 2.0 < y C < 4.5 and 1 < p C T < 20 GeV/c are measured for Υ(1S)D 0 and Υ(1S)D + final states at √ s = 7 and 8 TeV. The measured cross-sections are in agreement with DPS expectations and significantly exceed the expectations from the SPS mechanism. The differential kinematic distributions for ΥD 0 and ΥD + are studied and all are found to be in good agreement with the DPS expectations as the main production mechanism. The measured effective cross-section σ eff is in agreement with most previous measurements.