The role of multi-parton interactions in doubly-heavy hadron production

Beauty and charm quarks are ideal probes of pertubative Quantum Chromodymanics in proton–proton collisions, owing to their large masses. In this paper the role of multi-parton interactions in the production of doubly-heavy hadrons is studied using simulation samples generated with Pythia, a Monte Carlo event generator. Comparisons are made to the stand-alone generators BcVegPy and GenXicc. New methods of speeding up Pythia simulations for events containing heavy quarks are described, enabling the production of large samples with multiple heavy-quark pairs. We show that significantly higher production rates of doubly-heavy hadrons are predicted in models that allow heavy quarks originating from different parton–parton interactions (within the same hadron–hadron collision) to combine to form such hadrons. Quantitative predictions are sensitive to the modelling of colour reconnections. We suggest a set of experimental measurements capable of differentiating these additional contributions.


Introduction
The masses of the charm and beauty quarks are so high that their production at hadron colliders is dominated by perturbative QCD processes making them ideal probes of the partonic interactions that led to their formation. Hadrons containing two beauty and/or charm quarks are referred to as doubly-heavy hadrons and, here, include both mesons and baryons with a net heavy flavour such as B + c mesons and ++ cc baryons, and quarkonia with zero net flavour. Large samples of doubly-heavy hadrons have been collected at the Large Hadron Collider (LHC) and the study of their properties can provide a unique insight into the role of multi-parton interactions (MPIs) in hadron formaa e-mail: tom.hadavizadeh@cern.ch (corresponding author) tion. Because c and b quarks are too heavy to be produced non-perturbatively, doubly-heavy hadrons can only be formed by the effective coalescence of two perturbatively produced heavy quarks. This is fundamentally different from singly-heavy hadrons, in which the heavy quark is confined together with light quark(s) that can be produced non-perturbatively. The fact that heavy hadrons can be produced via MPI has been confirmed not only in inclusive measurements of heavy-hadron-pair cross sections [1][2][3][4] but also in differential measurements of J/ψ rates vs chargedtrack multiplicity [5]. However, these measurements do not by themselves show conclusively whether two heavy quarks from different parton-parton interactions can join to form hadrons. 1 Recent measurements of the newly-discovered doubly-charmed tetraquark T ++ cc indicate the production has similarities with processes involving different parton-parton interactions [6]. In this paper, we point out that this question can be addressed by considering non-onium doubly-heavy hadrons, such as B + c and ++ cc . Moreover, such measurements will place constraints on models of colour reconnections (CR), which are relevant to a broad range of hadroncollider physics studies [7][8][9][10].
Samples of doubly-heavy hadrons can be simulated using Monte Carlo event generators such as Pythia [11,12]. However, as these hadrons require the chance coalescence of two heavy particles into a bound state, their generation can be prohibitively slow. For this reason some rare doubly-heavy hadrons, for example, B + c , cc and bc , are at the moment usually generated using dedicated generators, such as BcVegPy [13] or GenXicc [14], that perform fixed-order matrix-element calculations. These generators are then interfaced with event generators to simulate the rest of the event evolution and hadronisation. However, (a) (b) (c) Fig. 1 Examples of production mechanisms for heavy quarks in proton-proton collisions. The incoming, outgoing and intermediate particles of the process considered to be the hardest process are highlighted in red. In the case of flavour excitation, theb quark shown at the bottom represents the companion quark produced as a result of the initial-state evolution the fixed-order matrix-element generators assume that both of the heavy quarks are produced in a single parton-parton interaction, for example gg → B + c bc or qq → B + c bc. This ignores the role that MPIs could play in generating heavy quarks that could contribute to the formation of such hadrons.
In this paper, we develop a method to enhance the rate of hadrons containing one or more heavy quarks in simulation samples generated with Pythia. This enables us to investigate the role of MPIs in the production of doubly-heavy hadrons in the context of inclusive event samples, in which the formation of such hadrons would normally be exceedingly rare. The generated samples are studied to identify properties that can differentiate between the contributions from single parton scattering (SPS) and double parton scatering (DPS) mechanisms. Simulations are carried out with Pythia version 8.306 using the default Monash tune [15] and at a proton-proton centre of mass energy of 13TeV. The hard interaction is simulated using Pythia 8.306, BcVegPy 2.2 or GenXicc 2.1 while the fragmentation process is simulated using the default Pythia Simple Shower framework [16,17]. Finally, measurements using LHC data that are able to shed light on the production mechanism are proposed and their feasibility discussed.
Throughout this paper charge conjugation is not implied when describing cross sections. The generator BcVegPy only produces B + c mesons, therefore to avoid ambiguity all cross sections are explicitly only produced for the explicitly specified charge. In general the cross section predictions from BcVegPy should be doubled when comparing to measurements of σ (pp → B ± c X ).

Sources of heavy quarks
In proton-proton collisions, the QCD production mechanisms for the heavy c and b quarks can be split into three categories referred to as pair creation, flavour excitation and gluon splitting in parton showers [18]. The processes are classified according to the interaction with the largest momentum transfer, referred to here as the hard interaction. Pair creation involves a gg → QQ 2 or qq → QQ hard interaction, as shown in Fig. 1a, that, in the absence of significant initial-state radiation, creates outgoing heavy quarks with equal and opposite transverse momenta. The resulting heavy hadrons formed from the heavy quarks similarly have a strong tendency to be back-to-back in the transverse plane, as shown in Fig. 2 for bb production.
Flavour excitation is the process involving one heavy quark: Qg → Qg or Qq → Qq, represented in Fig. 1b. In this process a virtual QQ pair is produced as part of the initial-state evolution of one of the incoming protons, and one of them, say the Q, subsequently interacts with a (non-heavy) parton from the other proton. TheQ (a.k.a. the "companion" quark of the scattered heavy quark [19]) is ejected as part of the initial-state evolution of the incoming remnant at a lower scale, with less transverse momentum and significantly less correlation with the direction of the Q, as shown in Fig. 2.
Heavy quarks can also be produced via gluon splittings during parton showers. A typical example would be a hard gg → gg interaction followed by a subsequent g → QQ splitting in the subsequent initial-or final-state shower evolution, as shown in Fig. 1. Although this figure shows one of the outgoing gluons from the hard interaction directly splitting to heavy quarks, that is just for simplicity; in principle any gluon produced within a shower above the heavy quarkmass threshold could result in heavy quarks. As gluon-gluon interactions have a large cross-section at the LHC, this constitutes a significant contribution to the heavy-quark production mechanisms. For final-state gluon splittings, the resulting QQ pair will be boosted in the direction of the parent gluon. Events in which two singly-heavy hadron are produced by this mechanism tend to have smaller angles between the two heavy hadrons, as shown in Fig. 2.

Sources of doubly-heavy hadrons
To create doubly-heavy hadrons that are not quarkonium states, two QQ pairs must be produced during the perturbative evolution of the collision. An example of an SPS mechanism contributing to this process is shown in Fig. 3a: hard bb pair creation followed by a g → cc splitting during the shower evolution. Equivalent processes involving flavour excitation or double gluon splitting within a single SPS are of course also possible.
When allowing for MPI, the two QQ pairs may also be produced in two different parton-parton interactions (still within the context of a single hadron-hadron collision). This is what we label DPS. Two examples, double pair creation and double flavour excitation, are shown in Fig. 3b and c respectively, again with other combinations of pair creation, flavour excitation, and/or gluon splittings obviously also possible. In these diagrams the two parton interactions have been highlighted in different colours to clarify the origin of the partons.
In events with more than two parton-parton interactions, SPS mechanisms could contribute from any one of the single parton-parton interactions, whilst DPS mechanisms could contribute from the combination of any two.
Once the appropriate quarks have been produced in the collision, only pairs that are sufficiently close in phase space and which have a non-zero probability to be in an overall colour-singlet state, have a chance to form an on-shell doubly-heavy hadron.
Measurements of the cross sections of multiple heavy hadrons suggest that MPIs play a significant role in the production of multiple heavy quark pairs at hadron colliders [1][2][3][4][5]. However, the question of how partons originating from different parts of the protons become bound into hadrons is still afflicted with significant uncertainties. In general-purpose event generators like Pythia, this is controlled by a combination of perturbative heavy-quark production mechanisms (hard scatterings, MPI, and parton showers) and semiempirical models of colour reconnections with [10,20] and without [7,8] space-time dependence. The simple diagrams in Fig. 3 demonstrate how B + c mesons formed from thebc combinations could provide an ideal probe into the hadronisation process. This is unique to doubly-heavy hadrons, since light quarks are mainly created nonperturbatively and hence do not have the same character of being associated with specific short-distance processes in the colliding protons.

Efficient simulation of events with heavy hadrons in Pythia
Generating unbiased events with multiple pairs of heavy quarks and doubly-heavy hadrons with Monte Carlo event generators can be very time consuming as few events will fulfil the requirements to form the doubly-heavy hadrons. A method of enhancing the efficiency to produce events containing heavy quarks in Pythia is outlined here, and can be applied to both singly-and doubly-heavy hadrons.
Pythia provides user-configurable classes called UserHooks aimed at allowing the user to inspect and veto events at different stages during the event evolution. These can be exploited to veto events that do not contain the requisite heavy quarks early on in the generation, removing time spent evolving and hadronising events that will never be accepted.
The UserHook stages that are utilised to improve the efficiency are: • Hard-process-level veto: This veto inspects the event after the most energetic parton interaction has occurred, as represented in Fig. 4a; • Event-evolution veto: In Pythia the event is evolved from the hard-interaction scale down to the hadronisation scale. During this process, the event can be inspected when the evolution reaches an arbitrary (user-defined) value of the evolution scale, illustrated in Fig. 4b; • Parton-level veto: Once the evolution has reached the hadronisation scale, the event can be inspected before any hadrons are created, leaving the complete partonic event as shown in Fig. 4c; • Post-hadronisation veto: It is possible to inspect the event after the hadronic states have been formed as shown in Fig. 4d.
Obviously, vetoes that can be placed early on can lead to far bigger efficiency gains than ones placed later, and a veto at the post-hadronisation level, after most of the event activity has already been treated (apart from decays and, optionally, hadronic rescatterings) will not make much of a difference at all.

Simulating final states involving one QQ pair
The creation of a heavy quark, Q, involves a physical momentum transfer of at least O(m Q ), regardless of whether the production mechanism is a hard process, MPI, or a g → QQ shower branching. Given an event-generator evolution algorithm that is ordered in a measure of momentum transfer, when the evolution reaches a scale of order m Q one could thus immediately veto any events that do not contain a specific desired number of heavy quarks. 3 Doing so effectively "saves" the time it would otherwise have taken to evolve such events from m Q to the perturbative cutoff, as well as the time required to hadronise them. Since the running value of α s is largest near the cutoff, where the number of evolving partons is also highest, and hadronisation can consume a further significant amount of time (in particular when colour reconnections are involved), a veto at the scale m Q should produce order-of-magnitude efficiency gains, compared with a standard simulation in which events are fully evolved and hadronised before the event is inspected. A caveat to implementing such a strategy for Pythia, however, is that its default transverse-momentum evolution scales [16,17] for heavy-quark processes are only guaranteed to be greater than m Q for hard processes and for ISR, whereas the probability densities for FSR g → QQ splittings and MPI gg → QQ processes are non-zero all the way down to p ⊥evol → 0. This is illustrated in Fig. 5, for a representative case of gg → gg hard-scattering events withp ⊥ = 25 ± 0.5 GeV in proton-proton collisions with √ s = 10 TeV. An interesting follow-up question for future work is thus whether it would be physically justifiable to reformulate these algorithms in terms of a measure that would associate all heavy-quark production processes with scales ≥ m Q . (This is, e.g., the choice made in the Vincia shower [21], which however still relies on Pythia's MPI evolution.) For now, we accept that this subtlety will force a tradeoff between realising the full possible efficiency gains and "missing" a small fraction of signal events, which we will comment on further below. We therefore phrase our implementation in terms of an arbitrary veto scale, not necessarily equal to the heavy-quark mass, which can be varied to determine if an "acceptable" tradeoff can be found, which may vary from application to application. Specifically, we define the following Userhook vetoes: Hard-process-level veto: At this stage, nothing is known about the subsequent shower or MPI evolution, except for what the starting scale for those evolutions will be. Our veto function only accepts events that fulfil at least one of the following two conditions: (1) the hard process itself contains the requisite heavy flavour (by which we include any onium containing it or a heavier quark that can decay to it), in which case a flag may also be set to bypass any downstream vetoes, or (2) the starting scale for MPI and showers is above our user-defined veto scale, so that we want to give MPI and/or showers a chance to produce the heavy flavour. This essentially means that gg → gg events withp ⊥ < O(m Q ) can be rejected already at this stage, with minimum processing.

Event-evolution veto:
If the hard-scattering process did not contain the requisite heavy flavour but was allowed a chance to produce it via MPI and/or showers, the event is inspected again when the evolution reaches our veto scale, and is now rejected if the required flavour (again including onia and/or heavier flavours) is still not present in the event.
The improvement in efficiency when generating samples with these two UserHooks is investigated for samples of events containing bb or cc. The time taken to generate the QQ pairs is compared to a baseline without the UserHooks included. All timing tests are performed using an Apple M1 MacBook Pro. 4 The relative speed-up and fraction of events missed due to the evolution scale definition are shown for bb pairs in Fig. 6. A significant improvement in efficiency is found when generating bb pairs with the UserHooks. The improvement is less significant when generating cc pairs because the smaller c-quark mass means the event evolution must continue further before the event can be vetoed.
The p T distribution of B hadrons in events that are not retained by the UserHooks are shown in Fig. 7. This sample, produced with the Simple Shower model misses bb pairs produced in both the parton shower and as additional MPI interactions. Overall, setting ap T scale of 4 GeV gives a factor 10 improvement in simulation speed, but leads to a small distortion in the p T spectra of the generated b hadrons.
The impact of these efficiency improvements can further contextualised in terms of the typical time taken to generate specific singly-heavy hadrons with Pythia. The typical times with and without the developed UserHooks are listed in Table 1.

Simulating final states involving multiple QQ pairs
When simulating events with multiple pairs of heavy quarks, the same Userhook vetoes described previously in Sect. 4.1 can be utilised. In principle, for B + c production and the like, it would be useful to apply an event-evolution veto first at O(m b ) and then again at O(m c ). However, with current versions of Pythia, the event can only be inspected at a single value of the evolution scale. Therefore, we set the the eventevolution threshold according to the heaviest quark being simulated, while the end-of-evolution parton-level veto is used to check for any required secondary heavy flavour.
Parton-level veto: If both a bb and a cc pair is requested (and/or onia containing them), the presence of the lighter of the two flavours is checked for at the end of the partonlevel evolution.
We believe that, by making minor modifications to the Pythia source code it could be possible to allow the eventevolution veto to inspect the event multiple times, at different scales during the event evolution (e.g. at both the b and c quark masses). This would save the time spent evolving the event from the lower mass scale to the hadronisation scale, providing benefits in addition to those demonstrated here.

Simulating final states with specific hadrons
Further UserHooks can be placed at the partonic event level to increase the efficiency of generating specific hadrons. For example, when generating B + c mesons, the invariant mass of  all c andb pairs (and vice versa) can be calculated to veto events in which there are no combinations that have a small enough invariant mass to form a hadron. The timing improvements when using the additional levels of UserHook are shown in Fig. 8, where the first UserHook to be included requires the smallest invariant mass of any bc orbc pair to be below 10 GeV/c 2 , the second requires the partonic-level event to contain both bb and cc pairs, and the third and fourth correspond to the hard process and event evolution vetoes already described, here applied just to the b quark. Using combinations of the UserHooks described above it is possible to improve the efficiency of generating both singly-heavy and doubly-heavy hadrons with Pythia. Typical times to generate doubly-heavy hadrons are listed in Table 2.
The UserHooks that inspect the event at the partonlevel or after do not result in any further missed hadrons, in contrast to the hard-process level and event-evolution vetoes discussed in Sect. 4.1. Therefore when generating doublyheavy hadron with these extra requirements, the fraction and kinematic distributions of the missed hadrons will be the same as for singly-heavy hadrons.

Comparison of doubly-heavy hadron kinematics
The efficiency improvements detailed in the previous section enable a range of comparisons that were previously unfeasible. For example, 250 000 B + c mesons have been generated for this study which would take approximately 10 CPU days with the UserHooks, or 140 CPU days without. In this section distributions of doubly-heavy hadrons generated with Pythia are compared to other standalone generators.

B + c mesons
The production of B + c mesons at hadron colliders is generally calculated in terms of the SPS process gg → B + c bc that is assumed to dominate [22][23][24], although there are some suggestions that flavour excitation processes could play an important role [25]. Currently, DPS production has not been considered. Measurements of the ratio of the B + c to the B + cross-section times branching fraction are found to be consistent with predictions of SPS contributions only [26], although the theoretical predictions have large uncertainties. In contrast, measurements of ϒ D cross section are consistent with the presence of DPS [4]. The ϒ D system also requires the creation of a bb and a cc pair, implying that there is the potential for B + c mesons to be similarly produced. Samples of B + c mesons are generated with Pythia and BcVegPy. The specific generation settings are listed in Table 3.
The kinematic distributions and differential cross-sections are compared in Fig. 9. The sample generated with Pythia is split according to the origin of the c-quark and b-quark that formed the B + c meson: those that originated in the same parton-parton interaction are categorised as SPS and those from different parton-parton interactions are categorised as DPS. This is determined by comparing the parent quarks to the record of parton interactions stored in Pythia's PartonSystems class. The production cross-section predicted by Pythia receives a significant contribution from DPS and is found to be larger than the production crosssection obtained from BcVegPy.
The kinematic distributions of the B + c alone do not provide significant discrimination power to the presence of hadron formation in DPS, other than in the absolute normalisation of the cross section. Measurements of the B + c cross section relative to B + c mesons have been performed [26][27][28], but extracting the absolute cross section requires theoretical predictions of the branching fractions.  A better way to investigate the role of SPS and DPS contributions is to measure the production cross section as a function of the number of parton-parton interactions in a collision. For hadrons formed in SPS processes, increasing the number of parton interactions would linearly increase the number of opportunities to form the hadron, as each new parton interaction would present one more opportunity for the hadron to form. However, hadrons formed in DPS processes would see the rate of formation increase quadratically with the number of interactions, as each hadron requires two parton interactions to form. These different relationships can be exploited to differentiate the components by considering the ratio of doubly-heavy to singly-heavy hadron cross sections, as a function of the number of parton-parton interactions. This ratio would be flat if singly-and doubly-heavy hadrons are produced by the same mechanism -SPS -while it would increase linearly if there is a nontrivial DPS component to doubly-heavy hadron production.
In Pythia, both mechanisms are present, while in BcVegPy, a single gg → B + c bc interaction is produced for each event which is then passed to Pythia for showering, MPI, and hadronisation. In this case, there is therefore no opportunity for heavy quarks from different parton-parton interactions to form the B + c meson and the production is independent of the total number of parton-parton interactions.
The cross-section ratio of B + c to B + mesons is compared for Pythia and BcVegPy in Fig. 10 as a function of the number of parton-parton interactions. In this figure no kinematic requirements have been placed on the rapidity or transverse momentum of the B + c meson or final-state particles. As expected, the contribution from DPS varies as a function of the number of parton interactions in the event. A significant enhancement is seen in events with many parton interactions. This quantity is not directly observable, but it is highly correlated to the number of particles in a collision. A more realistic, experimentally-accessible distribution is therefore the crosssection ratio as a function of the number of charged particles within a given acceptance. In the right pane of Fig. 10 and throughout we use the acceptance of the LHCb experiment (2.0 < η < 4.5), but other ranges would show similar results.
The distribution of the number of charged particles within 2.0 < η < 4.5 in events containing a B + c meson is shown in Fig. 11 for the different generation configurations. In the Pythia SPS and BcVegPy samples there are means of approximately 34 charged particles per event, whilst the Pythia DPS sample has a higher mean of around 41.
The significant difference in these distributions would allow measurements of the differential cross section with respect to B + mesons to differentiate between B + c mesons produced in SPS and DPS.

++ cc baryons
The doubly-charmed baryon ++ cc has been observed [29,30], and could similarly receive contributions from DPS production mechanisms. Samples of ++ cc baryons are simulated with both Pythia and GenXicc. In Pythia, the formation of baryons with multiple heavy quarks is only possible with the so-called QCD colour-reconnection option [8] which allows for the formation of doubly-heavy diquarks via string junctions [31]. The specific Pythia settings are listed in Table 4.
The GenXicc generator can produce different configurations of the initial charm hadrons. The processes gg → ++ cc cc and gc → ++ cc c are simulated. The kinematic distributions of the ++ cc samples are shown in Fig. 12, where again the Pythia samples have been split according to the whether the two heavy quarks originated from the same or different parton-parton interactions. Analogously to the case for B + c , a significantly larger cross section is predicted by Pythia, due to the presence of the DPS component.
The ratio of the differential cross-section between the ++ cc and D + hadrons are shown as a function of the number of parton interactions and number of charged particles within the pseudo-rapidity region 2.0 < η < 4.5 in Fig. 13. The crosssection ratio is clearly able to differentiate between DPS and SPS production mechanisms. The difference in the cross-  section ratios for the gg → ++ cc cc and gc → ++ cc c processes simulated by GenXicc may be caused by the different typical momentum transfers in the two processes. When simulating the rest of the underlying event with Pythia, this quantity determines the maximum scale of the event evolution, therefore the typically lower momentum transfers in gc → ++ cc c processes may cause fewer parton-parton interactions to be generated in the subsequent event evolution.

+ bc baryons
Similarly, samples of + bc can be produced with Pythia and GenXicc. These baryons have not yet been observed. The kinematic distributions of the + bc baryons in samples generated with Pythia and GenXicc are compared in Fig. 14, and the ratios of the differential cross section with respect to the B + meson are shown in Fig. 15 as a function of both the

Quarkonia
For the purposes of this study, quarkonia are considered doubly-heavy hadrons as they contain two heavy quarks. However, unlike the previously discussed doubly-heavy hadrons, quarkonia can be formed from a single heavy-quark pair. As a result, contributions to quarkonium production from heavy quarks in different parton interactions will only give a subleading contribution to the total production rate. This can be observed in Fig. 16 where the ratio of J/ψ to D + production cross sections are shown as a function of the number of parton interactions and the number of charged tracks within 2.0 < η < 4.5. There is a contribution from J/ψ mesons formed from c and c quarks from different parton interactions, however this is much smaller than the total rate of J/ψ production. To increase the fraction of events with DPS contributions, events with both a J/ψ meson and two additional charm hadrons can be reconstructed. This removes events in which there is only a single quark pair and leads to a different set of measurements that can be made as discussed in the next Section.

Associated production of singly-and doubly-heavy hadrons in events with multiple QQ pairs
In proton-proton collisions that produce two pairs of heavy quarks, i.e. cccc, ccbb or bbbb, information about the production mechanisms can be inferred from the relative properties of a doubly-heavy hadron and two singly-heavy hadrons that can be formed from the additional heavy quarks. Examples of the different combinations of doubly-and singly-heavy associated productions are listed in Table 5, where singlyheavy hadrons containing a heavy quark are represented by X Q . Only the B + c meson, b-and c-flavoured quarkonia and ++ cc baryon have currently been observed [29,30,[32][33][34][35].
Of the processes listed in Table 5, those categorised as doubly-heavy mesons and doubly-heavy baryons require quarks from different QQ pairs to combine in order to form. This configuration is referred to as mixed. In contrast, the quarkonium in ψ(nS)X b Xb or ϒ(nS)X c Xc events contain differently flavoured quarks compared to the associated

(b) (a)
singly-heavy hadrons, and therefore must be formed from an individual heavy-quark pair. This configuration is referred to as unmixed. In events containing a quarkonium, X Q and XQ hadrons where the flavours of the singly-heavy hadrons match the doubly-heavy flavour, the quarkonium could have been formed from a single QQ pair, i.e. unmixed, or from the combination of a Q andQ from two different QQ pairs, i.e. mixed. This is demonstrated for the charm quark in Fig. 17.
Therefore, for the B + c , cc and bc hadrons any indication that DPS processes contribute to the production would indicate that heavy quarks from different parton interactions can form hadrons. In contrast, the indications that DPS processes contribution to ψ(nS)X b Xb or ϒ(nS)X c Xc events implies that there are events with multiple heavy-quark pairs from different parton interactions, but doesn't give any information about the hadronisation of those different pairs. The quarkonia final states ψ(nS)X c Xc and ϒ(nS)X b Xb can give information about both the presence of multiple heavy quark pairs and hadronisation in DPS, but care must be taken to separate the mixed and unmixed contributions.
The kinematic distributions of the singly-and doublyheavy hadrons depend on the specific production mechanism and differ for those possible within SPS or DPS. Collisions in which both heavy hadron pairs originate from a SPS, for example the process shown in Fig. 3a, have correlations between the kinematics of at least one pair of heavy quarks. Processes that involve DPS introduce the possibility of double flavour excitation processes, such as Fig. 3c, in which the kinematics of the heavy quarks that form the hadron can be relatively uncorrelated to the remaining two companion heavy quarks. 5 The production of quarkonia is complicated further by contributions from colour-octet mechanisms [36][37][38], which are included by default in the Pythia simulation samples.

Studies with B + c X b Xc events
The relative kinematic distribution in B + c X b Xc events are studied to determine if there are differences between the DPS and SPS sub-samples for the mixed configuration, and to compare with the standalone generator BcVegPy. The angular separation in the transverse view between the B + c and the associated X b hadron φ(B + c , X b ) is plotted against the same quantity between the B + c and associated X c hadron φ(B + c , Xc) in Fig. 18 for the different generators and configurations. To ensure an unambiguous association between the heavy quarks, only events with a total of one B + c meson, one Xc hadron and one X b hadron are used.
In the sample generated with BcVegPy the generated events are found predominately with large angles between the B + c meson and X b hadron, and small angles between the B + c and Xc hadrons. This is consistent with the topology shown in Fig. 3a in which the bb pair are produced in the hardest interaction, and therefore back-to-back in the transverse plane. The c-quark resulting from a gluon splitting would then be produced in a cone around the B + c direction. In contrast, the DPS sample produced by Pythia introduces the possibility (c) (b) (a) Fig. 18 Angular separation in the transverse view between B + c mesons and the associated X b or Xc hadron in events generated with the BcVegPy and Pythia generators of additional production mechanisms including those shown in Fig. 3b and c. As such the distribution of events in the 2D plane is less localised as a result of the contributions from many different associated production correlations.

Studies with ϒ(1S)X c Xc events
Events containing both an ϒ(1S) meson and X c Xc pair can only receive contributions from the configuration referred to as unmixed. It therefore provides a suitable system to test whether MPIs contribute significantly to events with multiple pairs of heavy quarks, but cannot provide insight into the hadronisation of heavy quarks from different parton interactions.
The relative transverse distributions between the ϒ(1S) meson and the X c and Xc hadrons are shown in Fig. 19a and b. The two-dimensional distributions show a clear difference in the relative distributions of the two types of process. In the SPS process the ϒ(1S) meson has a strong tendency to be produced back-to-back to both charm hadrons in the transverse plane. This could result from gg → ϒ(1S)g parton interactions in which the outgoing gluon subsequently splits g → cc. In contrast, in DPS processes there is little correlation between the transverse directions. The predictions can be compared to measurements by LHCb [4] by making a one-dimensional projection, as shown in Fig. 19c. Only events in which both the ϒ(1S) and one of the corresponding X c hadron are within the LHCb acceptance are compared in this figure. The data are consistent with the predictions including DPS, but the significant difference between the SPS and DPS samples is diluted when projected onto one dimension. This strongly motivates performing measurements in which both associated hadrons are reconstructed, such that the two dimensional distributions can be determined. Finally, the ratio of the differential cross section relative to the B + meson is shown in Fig. 19d as a function of the number of charged particles within 2.0 < η < 4.5. It similarly shows a strong separation power between the two process types.

Studies with J/ψ X c Xc events
Samples containing a J/ψ meson and two charm hadrons are generated with Pythia. The samples are split according to whether the J/ψ meson was created in a single interaction, or from quarks in different parton interactions, as before. Additionally, the ancestors of the charm quarks are studied to determine if the two quark lines are mixed or unmixed, as defined in Fig. 17.
The ratio of the differential cross sections with respect to the D + meson is compared for the four categories, as shown in Fig. 20. The two DPS categories show an increasing trend with both the number of parton interactions and charged particles within 2.0 < η < 4.5. The fraction of events produced in DPS processes is significantly increased with respect to the samples in which just the J/ψ was reconstructed in Sect. 5.4.
The relative transverse directions of the three particles is shown for the four categories in Fig. 21. The two unmixed categories show similar features to those already seen in ϒ(1S)X c X c events. In contrast the two mixed categories have distinct distributions, with the DPS mixed category having a slight tendency to have the two charm quarks produced in the same direction as the J/ψ , a topology that may arise due to events with two flavour excitation interactions.
The presence of these events categorised as DPS mixed would indicate that heavy quarks from different parton interactions can hadronise to form a J/ψ meson. Using a combination information from the relative cross section ratio and relative transverse directions, the presence of this contribution could be determined.

Studies with ϒ(1S)X b Xb events
Similar studies can be performed on the ϒ(1S)X b X b system. Large samples are generated with Pythia requiring a minimal event content of bbbb. Events with one ϒ(1S) meson and a X b X b pair are selected, and the relative distributions studied. The events are split according to whether the One-dimensional projections of the relative angular distributions, compared to measurements from Ref. [4]. (Bottom right) Relative differential cross-section of ϒ(1S)X c Xc events with respect to B + as a function of the number of charged particles within the pseudo-rapidity region 2.0 < η < 4.5 Fig. 20 Ratio of differential cross-section of J/ψ X c X c and D + hadrons as a function of (left) the number of parton interactions in a collision and (right) the number of charged particles within the pseudo-rapidity region 2.0 < η < 4.5, as generated with Pythia b and b that formed the ϒ(1S) were from the same or different parton systems. As before, the ancestry of the b and b quarks is studied to determine if the configuration is mixed or unmixed. The relative differential cross-section relative to the B + meson is shown in Fig. 22 as a function of the number of parton interactions and the number of charged particles within 2.0 < η < 4.5. A similar trend is observed in this system.
Additionally, the relative transverse directions of the ϒ(1S), X b and X b hadrons is studied for the four production categories. The distributions, shown in Fig. 23 show a similar distribution for the SPS unmixed as before. In this system, events categorised as DPS unmixed have a slight ten-dency to have an anti-correlation between φ(ϒ(1S), X b ) and φ(ϒ(1S), Xb). This could correspond to events in which one parton interaction produces the ϒ(1S) meson and another produces a bb pair from pair production that are back-to-back in the transverse plane. Therefore as φ(ϒ(1S), X b ) increases, the corresponding value of φ(ϒ(1S), Xb) decreases. The DPS unmixed distribution has a distinct distribution, potentially allowing discrimination.
6.5 Impact of the colour reconnection model As mentioned briefly in Sect. 5.2 and discussed in more detail in [8,9,12], the QCD CR model of [8] allows for coher-  [31]. In the limit of a small invariant mass between two heavy quarks, the junction topology (when allowed by the CR selection rules [8]) reduces to a doubly-heavy diquark. In Pythia, this mechanism is essentially the sole means by which doubly-heavy baryons can be produced at all, making the rates and spectra of such baryons -and their dependence on event characteristics such as N Charged -particularly sensitive probes of this type of colour-space ambiguities. (Accord-ingly, we note that all of our Pythia results for doubly-heavy baryons in this work were produced with the QCD CR option; otherwise the rates would be essentially zero.) The choice of CR model also has an impact on the relative size of the DPS contribution to doubly-heavy meson production. The effect of using the QCD CR model described in Table 4 on the slope of the ratio of B + c to B + differential cross sections is shown in Fig. 24. With the QCD CR model, the DPS contribution gets smaller, altering the slope of the Similarly, the relative cross section distributions are found to vary in J/ψ X c Xc samples generated with the QCD CR options enabled, as shown in Fig. 25. The contribution from the DPS mixed configuration decreases whilst DPS unmixed increases, implying that the likelihood for the heavy quarks from different parton-parton interactions to be combined into a single hadron is sensitive to the choice of CR scheme.
Additionally, the relative transverse direction distributions are found to differ for the DPS configurations. The corresponding distributions are shown in Fig. 26. Measurements of J ψ X c Xc events may additionally help to differentiate between the different models of colour reconnection.

Experimental measurements and feasibility
In the Pythia simulation studies performed for this paper, the production of doubly-heavy hadrons is predicted to have a significant contribution from DPS production processes. New measurements of the relative cross section for the doubly-heavy hadrons with respect to singly-heavy hadron as a function of the collision multiplicity would help identify if such contributions are present in nature, as proposed in Sect. 5. Unlike recent observations of strangeness enhancements in the ratio of B 0 s to B 0 cross sections [39], the enhancements from DPS are not expected to be localised. The most suitable doubly-heavy hadron for this would be the B + c meson. The significant yields reported in a selection of different papers are listed in Table 6. Studies may also be feasible for ++ cc baryons. The predicted fraction, , of B + c mesons produced in DPS processes varies as a function of p T (Fig. 27), implying the effects would be most pronounced at lowp T . This would motivate measuring the relative cross sections as a function of the number of tracks in different p T regions. The DPS fraction is not found to vary as a function of rapidity.
The contributions from DPS production mechanisms can also be studied in events with one quarkonium and two singlyheavy hadrons, as discussed in Sect. 6. These final states have the advantage that quarkonia can be efficiently reconstructed using leptonic final states. The singly heavy hadrons could be reconstructed exclusively, using a range of different final states, or inclusively, by exploiting the typical topologies of heavy-meson decays.  Table 4 (triangles) and with the default CR options (circles) Fig. 25 Ratio of differential cross-section of J/ψ X c X c and D + hadrons as a function of (left) the number of parton interactions in a collision and (right) the number of charged particles within the pseudo-rapidity region 2.0 < η < 4.5, as generated with Pythia both with the alternative CR options specified in Table 4 (triangles) and with the default CR options (circles) Fig. 26 Relative transverse distributions in J/ψ X c Xc events for DPS processes in both the mixed and unmixed configurations, as generated with Pythia with the colour reconnection options specified in Table 4 (b) (a)

Exclusive reconstruction of associated singly-heavy hadrons
For singly-heavy hadrons containing c-quarks the branching fractions of the main experimentally-efficient decay channels constitute a reasonable fraction of the total decay width, as listed in Table 7. Therefore it may be feasible to exclusively reconstruct a sufficient fraction of the charm hadrons to make the proposed measurements. Indeed, the cross-section of LHCb ++ cc → + c K − π + π + 1598 Run2 [43] LHCb ++ cc → + c π + 616 Run2 [43] Fig. 27 Fraction of B + c decays predicted to be produced by DPS processes as a function of (left) p T and (right) rapidity in simulations samples produced by Pythia Table 7 Branching fractions of experimentally efficient decay channels for different charm hadron species from Ref. [44] Decay mode B(%) D 0 → K − π + 3.95 ± 0.03 D 0 → K − π + π − π + 8.22 ± 0.14 D + → K − π + π + 9.36 ± 0.16 D + s → K + K − π + 5.39 ± 0.15 D + s → π + π − π + 1.08 ± 0.04 + c → pK − π + 6.28 ± 0.32 quarkonia plus one singly-heavy charm hadron have already been performed [1,4]. However, for hadrons contain b-quarks the typical decay channels have much smaller branching fractions and may decay further to charm hadrons, therefore the fraction of events containing a quarkonia in which it is possible to exclusively reconstruct two X b hadrons is small.

Inclusive reconstruction of associated singly-heavy hadrons
It is possible to inclusively reconstruct singly-heavy hadrons by taking advantage of the common features of many decay channels, namely a secondary decay vertex that is significantly displaced from the primary interaction. In order to measure the kinematic relationships between the singly-and doubly-heavy hadrons the position of displaced vertices must be measured, along with a flavour tag labelling the vertex as a b-or c-hadron. Alternatively, the singly-heavy hadrons could be reconstructed as charm and beauty jets. In this case the reconstruction of the full jet may allow a better approximation of the parton-level kinematics to be determined rather than the heavy hadron kinematics. Jet flavour-tagging algorithms use similar displaced vertex tagging algorithms as already discussed to determine if a jet is the result of a heavy flavour quark. Recent measurements of heavy dijets at LHCb claim a reconstruction efficiency for bb and cc dijets of about 15% and 1.5%, respectively [45]. Therefore in a significant fraction of events containing quarkonia it may be feasible to additionally reconstruct these jets, assuming that the rate of fake jets can be sufficiently controlled.

Conclusions
Studies have been performed using the Pythia Monte Carlo event generator to investigate the contribution from DPS to the production of doubly-heavy hadrons. New UserHooks have been developed for Pythia to produce events with one or more heavy quarks more efficiently, enabling the samples to be produced in a realistic time frame. Comparisons have been made to SPS predictions from BcVegPy and GenXicc. The studies show significant and measurable differences in the ratio of cross sections between doubly-and singly-heavy hadrons as a function of the collision multiplicity that can discriminate between SPS and DPS production mechanisms. Further, the DPS mechanism is sensitive to the assumed scenario for colour reconnections. A range of measurements have been proposed, including cross-section ratios for B + c and ++ cc decays as well as a combination of cross-section ratios and relative transverse directions in events containing a quarkonium and two singly-heavy hadrons. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.