1 Introduction

The study of angular correlations of heavy-flavour (charm and beauty) particles in ultra-relativistic hadronic collisions allows the investigation of fundamental properties of quantum chromodynamics (QCD) in the heavy-flavour domain [1, 2]. In particular, the angular-correlation function between prompt D mesons and charged particles in proton–proton (pp) collisions is sensitive to the mechanisms of charm-quark production, fragmentation, and hadronisation into charm hadrons. The term “prompt” refers to D mesons originating both from direct charm-quark fragmentation and from the strong decay of excited charm resonances, and excludes D mesons produced from beauty-hadron weak decays. The typical structure of the two-dimensional correlation function between “trigger” D mesons and “associated” charged particles, expressed in terms of the pseudorapidity difference (\(\Delta \eta = \eta _{\mathrm{ch}}-\eta _{\mathrm{D}}\)) and azimuthal angle difference (\(\Delta \varphi = \varphi _{\mathrm{ch}}-\varphi _{\mathrm{D}}\)), features a near-side (NS) peak, centred at (\(\Delta \varphi ,\Delta \eta ) = (0,0)\), and an away-side (AS) peak at \(\Delta \varphi = \pi \) that is elongated along \(\Delta \eta \) [3]. Both peaks sit on top of an approximately flat continuum extending over the full (\(\Delta \varphi ,\Delta \eta \)) range. The height and width of the two correlation peaks are sensitive to the charm-quark production mechanisms.

Due to their large mass, the production of charm-quark pairs occurs through hard parton–parton scattering processes with large momentum transfers, and can be described by perturbative QCD (pQCD) calculations. While in leading-order (LO) processes the two quarks are produced back-to-back in azimuth, the next-to-leading-order (NLO) production mechanisms, such as flavour excitation and gluon splitting, can break this topology and alter the shape of the two correlation peaks [4]. Recent studies at the LHC suggest a relevant contribution from gluon splitting to heavy-quark production, possibly underestimated by Monte Carlo (MC) event generators with LO or NLO accuracy [5, 6]. The analysis of the properties of the near-side peak also allows for detailing the fragmentation and hadronisation processes which, starting from a coloured charm quark, lead to the formation of a D meson surrounded by a spray of colourless particles, experimentally identifiable as a charm jet. The production cross section of jets containing D mesons, as well as the jet momentum fraction carried by the D meson along the jet-axis direction, were recently measured by the ALICE and ATLAS Collaborations [7, 8]. In this regard, a systematic and differential analysis of the near-side correlation peak in terms of transverse momenta of trigger D meson (\(p_{\mathrm{T}}^{\mathrm{D}} \)) and other associated fragmenting particles (\(p_{\mathrm{T}}^{\mathrm{assoc}} \)), and of the angular distance of associated particles from the D mesons, can provide additional information with respect to measurements that treat charm jets as a whole entity.

In recent years, the study of high-multiplicity pp collisions has become of particular interest. The ALICE Collaboration has measured a faster-than-linear increase of prompt D-meson self-normalised yields for increasing relative event multiplicity in pp collisions at a centre-of-mass energy of \(\sqrt{s} = 7\) TeV, employing both central- and forward-rapidity multiplicity estimators [9]. A similar behaviour was also seen for open-beauty and hidden-charm hadrons, pointing towards sensitivity of heavy-quark production processes to event multiplicity. Complementary information on a possible dependence of charm-quark fragmentation and hadronisation on event multiplicity can be provided by the study of D-meson and charged-particle azimuthal correlations as a function of the event multiplicity, searching in particular for modifications of the near-side peak structure. This measurement is also crucial to validate the assumptions adopted to measure the elliptic-flow coefficient of D mesons and heavy-flavour decay muons in high-multiplicity pp collisions at \(\sqrt{s} = 13\) TeV by CMS [10] and ATLAS [11], respectively. In these measurements, the elliptic-flow coefficient is extracted from two-particle correlation function of such heavy-flavour particles with charged particles, and the short-range correlation peaks related to heavy-quark fragmentation are removed from the correlation function exploiting low-multiplicity events, assuming independence of the correlation peaks from the event multiplicity.

MC event generators, like PYTHIA [12, 13], EPOS [14, 15], or HERWIG [16,17,18], or pQCD calculations coupled to parton-shower software, such as POWHEG [19, 20], are widely used to reproduce ultra-relativistic hadronic collisions and provide predictions for a wide variety of physics observables. As discussed in Ref. [21], depending on the treatment of the various collision stages and implementation of specific features in each generator, such as hard parton–parton scattering matrix elements, parton-showering model, hadronisation algorithm, and underlying event generation, different predictions for D-meson and charged-particle correlation function will be obtained. A comparison of the predicted features of the correlation observables, in particular the peak yields and widths, with data measurements, can allow for validating and setting constraints to the MC generators [22].

A proper understanding of heavy-flavour correlations in pp collisions is also crucial in view of future studies in ultra-relativistic heavy-ion collisions. In the first stages of such collisions a deconfined state of strongly-coupled matter, the quark–gluon plasma (QGP), is created. While traversing the medium, heavy quarks interact with the QGP constituents through radiative and collisional processes [23, 24], losing energy and having their original directions modified. This is expected to lead to a modification of the angular-correlation function between final-state heavy-flavour hadrons and associated charged particles, with respect to that observed in pp collisions. Quantifying such modifications allows for investigation of specific properties of the QGP and its dynamics [2]. In particular, the angular-correlation function is sensitive to the relative contributions of the two energy-loss processes, and can shed light on the path-length dependence of energy loss [25,26,27,28]. Some first indications in this direction were provided by a recent measurement of \(\mathrm D^{0} \)-meson and charged-hadron angular correlations in gold–gold collisions at a centre-of-mass energy per nucleon pair of \(\sqrt{s_{\mathrm {NN}}} = 200\) GeV performed by the STAR Collaboration at RHIC [29]. Validating MC event generators against the correlation function of heavy-flavour particles measured in pp collisions is thus fundamental for a correct understanding of the same observables that will be measured in the future in lead–lead (Pb–Pb) collisions at the LHC.

In this article, ALICE measurements of azimuthal correlations of prompt \(\mathrm D^{0} \), \(\mathrm D^{+} \), and \(\mathrm D^{*+} \) mesons, together with their charge conjugates, with associated charged particles in pp collisions at \(\sqrt{s} \!=\! 13\) TeV at midrapidity are reported. For prompt \(\mathrm D^{0} \)-meson triggers the results of the correlation analysis are also reported as a function of the charged-particle event multiplicity, measured at forward and backward rapidity. With respect to previous ALICE publications in pp collisions [3, 21], the new measurements extend the \(p_{\mathrm{T}} \) range of D mesons up to 36 GeV/\(c\), and significantly improve the precision of the measured observables in the common \(p_{\mathrm{T}} \) range. Additionally, the measurements presented here, along with previous ALICE results at \(\sqrt{s} \!=\! 5.02\) and 7 TeV, enable the study of the possible evolution of the correlation distributions and of the peak features as a function of the collision centre-of-mass energy.

The article is organised based on the following structure. Section 2 describes the ALICE apparatus, as well as the data and MC samples used for this study. Section 3 highlights the procedure followed to obtain the azimuthal correlation function and to extract physical observables from it. In Sect. 4, the sources of systematic uncertainties affecting the results are detailed. In Sect. 5, the analysis results are presented and discussed, and a comparison with various model predictions is shown. Further model studies highlighting the specific contributions to the correlation function from initial- and final-state radiation and multi-parton interactions are reported in Sect. 6. A summary of the paper and its physics message is outlined in Sect. 7.

Table 1 The percentiles of the INEL>0 cross section of the four V0M-based event-multiplicity classes and the corresponding midrapidity charged-particle multiplicities. Systematic uncertainties on the charged-particle multiplicity values, derived from [38], are also reported
Full size table

2 ALICE detector, data, and MC samples

A complete description of the ALICE detector and its performance can be found in Refs. [30, 31]. The reconstruction of D mesons and charged particles was performed using detectors installed in the central barrel, with a pseudorapidity coverage of \(|\eta | < 0.9\) and a magnetic field of \(B = 0.5\) T parallel to the beam axis. In particular, the Inner Tracking System (ITS) [32] and the Time Projection Chamber (TPC) [33] were employed for the reconstruction of charged tracks and of primary and secondary vertices. The TPC, together with the Time-of-Flight (TOF) detector [34], also provided charged-particle identification (PID) information. The analysis also relied on detectors located along the beam line, at forward and backward rapidity. The V0 detector [35] is a set of scintillators covering the pseudorapidity ranges \(2.8< \eta < 5.1\) (V0A) and \(-3.7< \eta < -1.7\) (V0C), used for triggering, background-event rejection, and event-multiplicity estimation. The T0 detector is an array of Cherenkov counters, located along the beam line, at a distance of \(+370\) cm (T0A) and \(-70\) cm (T0C) from the nominal interaction point, and provides the collision starting time needed by the TOF [36].

The data sample used for the analysis consisted of pp collisions at \(\sqrt{s} = 13\) TeV collected during the 2016, 2017, and 2018 data taking periods, with a total integrated luminosity of about 29 nb\(^{-1}\), based on the visible cross section measured with the V0 detector [37]. The collisions were recorded if they satisfied a minimum-bias trigger, requiring the presence of signals in both V0 detectors in coincidence with a bunch crossing in the ALICE interaction region. This trigger was fully efficient for selecting events containing D mesons with \(p_{\mathrm{T}} > 1\) GeV/\(c\). Contamination of tracks from pile-up events (multiple collisions occurring in the same bunch crossing) was suppressed by discarding events where multiple primary vertices were reconstructed with the Silicon Pixel Detector (SPD), which constitutes the first two layers of the ITS. Timing information provided by the V0, as well as the correlation between the number of hits and the number of track segments in the SPD, were employed to reject beam–gas interactions. Only events with a reconstructed primary vertex within \(\pm 10\) cm from the nominal centre of the ALICE detector along the beam direction were considered to grant a uniform acceptance for the central-barrel detectors.

The multiplicity-differential analysis was performed in four independent multiplicity classes, defined in terms of the total energy deposit in the V0 detectors by charged particles passing through them (V0M amplitude). The rapidity gap between the V0 detectors and the central barrel in which the D mesons and charged particles were reconstructed assured the absence of significant auto-correlations between the correlation peaks and the multiplicity estimate. The V0M amplitudes were converted to percentiles of the inelastic collisions with at least one charged particle produced in \(|\eta | < 1\) (INEL>0), corresponding to about 75% of the total inelastic cross section, as described in Ref. [38].

The corresponding INEL>0 percentiles (\(\sigma \)/\(\sigma _{\mathrm{INEL>0}}\)) of the four V0M multiplicity classes are reported in Table 1, together with the related average number of charged particles, \(\langle \mathrm{d}N_{\mathrm{ch}}/\mathrm{d}\eta \rangle \) in \(|\eta | < 0.5\). A specific high-multiplicity trigger (V0HM) was used for the V0M multiplicity class I, to enhance the statistical precision of this particular class. The V0HM trigger recorded only events with a multiplicity large enough to pass a threshold of V0M amplitude. This trigger covered the whole span of V0M multiplicity class I. Only the data periods granting a uniform efficiency of the V0HM trigger inside the range covered by V0M multiplicity class I were considered for the multiplicity-differential analysis, resulting in an integrated luminosity specific to the V0HM trigger of about 7.7 pb\(^{-1}\).

To evaluate the corrections to the azimuthal-correlation measurements, several MC simulations of pp collisions at \(\sqrt{s} = 13\) TeV were used, produced with the PYTHIA 6.4.25 event generator [12] with the Perugia-2011 tune [39]. For the corrections specific to D mesons, a sample of pp collisions was produced with the same generator, with each event containing either a \(\mathrm{c}\overline{\mathrm{c}}\) or a \(\mathrm{b}\overline{\mathrm{b}}\) pair in the rapidity range \([-1.5,1.5]\). In addition, simulated events satisfying a minimum threshold of midrapidity charged-particle multiplicity were employed, as they provided sufficient statistical precision to evaluate the corrections for the V0M multiplicity class I. The simulations included the full description of the detector geometry, response, and conditions during the data taking via the GEANT3 package [40].

3 Analysis overview

The procedure for the evaluation of the D-meson (intending, in the context of this study, \(\mathrm D^{0} \), \(\mathrm D^{*}(2010)^{+} \), and \(\mathrm D^{+} \) mesons) and charged-particle azimuthal correlation function and the related corrections is described in Sect. 3.1 for the multiplicity-integrated analysis. The multiplicity-differential analysis largely followed the same approach, although some of the quantities and the corrections were evaluated independently in each V0M multiplicity class, or with a slightly modified procedure. Such differences are highlighted in Sect. 3.2. The extraction of physical observables from the fully-corrected correlation function, in common between the two studies, is discussed in Sect. 3.3.

3.1 Evaluation and correction of the azimuthal correlation function for the multiplicity-integrated analysis

All stages of the analysis were mostly unaltered with respect to those performed for the same study in pp collisions at \(\sqrt{s} = 5.02\) TeV, and are comprehensively described in Ref. [21]. Thus, they will be only briefly summarised in the following.

\(\mathrm D^{0} \)-, \(\mathrm D^{*}(2010)^{+} \)-, and \(\mathrm D^{+} \)-meson candidates and their charge conjugates, used as trigger particles in the analysis, were reconstructed from their hadronic decay channels \(\mathrm{D^{0} \rightarrow \mathrm{K}^{-}\pi ^{+}}\), with branching ratio BR = (3.95 ± 0.03)%, \(\mathrm{D^{+} \rightarrow \mathrm{K}^{-}\pi ^{+}\pi ^{+}}\) with BR = (9.38 ± 0.16)%, and \(\mathrm{D^{*+} \rightarrow \mathrm D^{0} \pi ^{+}} \rightarrow \mathrm{K}^{-}\pi ^{+}\pi ^{+}\) with BR = (2.67 ± 0.03)% [41]. A topological selection, exploiting the characteristic displacement of the D-meson decay vertices with respect to the primary vertex, and particle-identification information on the D-meson decay products were employed to suppress the combinatorial background. Further details on the selection are provided in Ref. [42]. The same criteria were followed in this analysis, apart from an optimisation of the selection values performed on the specific samples used to further increase the signal-to-background ratio of D-meson candidates. A fit to the invariant-mass distribution of selected D-meson candidates was performed as described in Ref. [42], in order to extract the D-meson yield, \(S_{\mathrm {peak\ region}}\), in a \(\pm 2 \sigma \) region from the centre of the invariant-mass signal peak, where \(\sigma \) is the width of the Gaussian component of the fit function describing the signal peak. The associated particles included charged primary [43] pions, kaons, protons, electrons, and muons with \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\) and \(|\eta | < 0.8\). The decay products of the trigger D meson were excluded from the associated-particle sample. The tracks reconstructed in the ALICE central barrel were accepted as associated particles if they satisfied selection criteria based on the quality of their reconstruction in the ITS and TPC detectors, as detailed in Ref. [3]. Additionally, a maximum distance of closest approach (DCA) of the track to the primary vertex of 1 cm in the transverse (xy) plane and along the beam line (z-direction) was required. This selection suppressed the contamination of non-primary particles to about 5% for \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), falling to below 1% for \(p_{\mathrm{T}}^{\mathrm{assoc}} > 2\) GeV/\(c\). As a result of the applied selection criteria, the associated-track reconstruction and selection efficiency ranged from about \(80\%\) for \(p_{\mathrm{T}}^{\mathrm{assoc}} = 0.3\) GeV/\(c\) up to about \(90\%\), increasing with \(p_{\mathrm{T}}^{\mathrm{assoc}} \).

Selected D-meson candidates within \(\pm 2 \sigma \) from the centre of the invariant-mass signal peak (“peak region”) were correlated with associated particles reconstructed and selected in the same event. A two-dimensional angular-correlation function, \(C(\Delta \varphi , \Delta \eta )_{\mathrm {peak\ region}}\), was built for each of the five D-meson \(p_{\mathrm{T}} \) intervals, ranging from 3 to 36 GeV/\(c\), and for the associated track \(p_{\mathrm{T}} \) interval \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\) and its sub-ranges: \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), \(1<p_{\mathrm{T}}^{\mathrm{assoc}} < 2\) GeV/\(c\), and \(2<p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\). The limited detector acceptance and efficiency for the reconstruction and selection of D-meson candidates and associated particles were accounted for by weighting each correlation pair by 1/\((A\times \epsilon )^{\mathrm{assoc}} \times 1\)/(\(A\times \epsilon )^{\mathrm{trig}}\), where A and \(\epsilon \) represent the acceptance and efficiency factors, respectively, evaluated using MC simulations. The \((A \times \epsilon )^{\mathrm{trig}}\) values were dependent on the event multiplicity. In particular, the D-meson selection efficiency decreased for low-multiplicity events, due to the degraded resolution on the primary-vertex position, which enters into several topological selection criteria. To account for this dependency, the \((A \times \epsilon )^{\mathrm{trig}}\) weights were evaluated and applied in narrow intervals of SPD tracklet multiplicity. The entries of the invariant-mass distributions of D-meson candidates were also scaled by 1/(\(A\times \epsilon )^{\mathrm{trig}}\) to allow a correct per-trigger normalisation of the correlation function, as detailed later. Additional losses due to pair acceptance effects were taken into account by applying a mixed-event correction. Specifically, D-meson candidate triggers were correlated with associated charged particles from other events with similar midrapidity event multiplicity and primary vertex location along the beam axis. Track segments reconstructed by associating hits in the two SPD layers and pointing to the reconstructed primary vertex (SPD tracklets) were used as the midrapidity multiplicity estimator for the event classification. In this way, a mixed-event correlation function, \(\mathrm {ME}(\Delta \varphi , \Delta \eta )_{\mathrm {peak\ region}}\), was built and used to weight the same-event correlation function \(C(\Delta \varphi , \Delta \eta )_{\mathrm {peak\ region}}\).

The correlation function \(C(\Delta \varphi , \Delta \eta )_{\mathrm {peak\ region}}\) also included a contribution from background D-meson candidates. This contribution was statistically removed by employing a sideband-subtraction technique. A sideband-region correlation distribution, \(C(\Delta \varphi ,\Delta \eta )_\mathrm{{sidebands}}\), was built by considering as trigger particles \(\mathrm D^{0} \)- and \(\mathrm D^{+} \)-meson candidates 4–8\(\sigma \) from the centre of the invariant-mass signal peak, in both directions, and \(\mathrm D^{*+} \)-meson candidates 5–10\(\sigma \) to the right of the invariant-mass signal peak centre. A mixed-event correction was applied to the sideband-region correlation function, following the same procedure described above for \(C(\Delta \varphi , \Delta \eta )_{\mathrm {peak\ region}}\). Subsequently, the sideband-region correlation function was subtracted from that of the peak region to obtain the signal correlation function, \(C(\Delta \varphi , \Delta \eta )_{\mathrm {signal}}\). The above procedure is described in Eq. 1

$$\begin{aligned}&\tilde{C}_\mathrm{{signal}}(\Delta \varphi , \Delta \eta ) = \frac{1}{S_\mathrm {peak\ region}}\left( \frac{C(\Delta \varphi , \Delta \eta )}{\mathrm {ME}(\Delta \varphi , \Delta \eta )} \bigg |_\mathrm {peak\ region} \right. \nonumber \\&\quad \left. - \frac{B_\mathrm {peak\ region}}{B_\mathrm {sidebands}}\frac{C(\Delta \varphi , \Delta \eta )}{\mathrm {ME}(\Delta \varphi , \Delta \eta )} \bigg |_\mathrm {sidebands}\right) , \end{aligned}$$
(1)

where the factor \(1/S_{\mathrm {peak\ region}}\) provides a per-trigger normalisation to the signal correlation function, as denoted by a \(\tilde{C}\) symbol. The terms \(B_{\mathrm {peak\ region}}\) and \(B_{\mathrm{sidebands}}\) quantify the number of background D-meson candidate triggers in the invariant-mass peak region and sidebands, respectively.

The two-dimensional correlation function \(\tilde{C}_{\mathrm{{signal}}}(\Delta \varphi \!,\!\Delta \eta )\) was integrated in the range \(|\Delta \eta | < 1\), obtaining the per-trigger azimuthal correlation function \(\tilde{C}_\mathrm{{signal}}(\Delta \varphi )\), in order to grant sufficient statistical precision. For the same reason, due to its symmetry around \(\Delta \varphi = 0\) and \(\Delta \varphi = \pi \), the correlation function was restricted to the \(0 \le \Delta \varphi \le \pi \) interval, averaging the symmetric points in the ranges [0, \(\pi \)] and [\(-\pi \), 0].

The residual contamination of non-primary associated tracks satisfying the DCA selection criteria was removed by multiplying \(\tilde{C}_\mathrm{{signal}}(\Delta \varphi )\) by a \(\Delta \varphi \)-differential purity correction factor \(p_{\mathrm{prim}}(\Delta \varphi )\), evaluated from MC simulations. A slight increase with \(p_{\mathrm{T}}^{\mathrm{assoc}} \) was observed for the \(\Delta \varphi \)-averaged value of this factor, which ranged between 0.94 and 0.99. Modulations as a function of \(\Delta \varphi \) up to 2% were obtained. The fraction of prompt D mesons (\(f_{\mathrm{prompt}}\)) accounted for approximately \(90\%\) of D mesons accepted by the topological and PID selection, with a slight increase for increasing \(p_{\mathrm{T}} \). The remaining contribution was composed of D mesons produced by beauty-hadron decays (feed-down D mesons). Thus, the azimuthal correlation function \(\tilde{C}_{\mathrm{{signal}}}(\Delta \varphi )\) included a contribution from feed-down D-meson candidate triggers. For small \(\Delta \varphi \) values, the shape of this contribution was distorted by the topological selection applied to the D-meson candidates, which was more efficient in selecting beauty-hadron decays featuring a small opening angle between the D-meson candidate trigger and the other decay particles. The natural shape of the feed-down contribution to the azimuthal correlation function was recovered by evaluating the amount of the distortion via MC studies, and applying a correction factor \(b_{\mathrm{B}-\mathrm{bias}}(\Delta \varphi )\) to the data correlation function, as explained in detail in Ref. [21]. The correction amounted to a maximum of 4.5% for \(\Delta \varphi = 0\) and was substantially smaller for larger \(\Delta \varphi \) values. After applying this correction, the feed-down contribution to the measured correlation function was removed as follows. A template of the per-trigger azimuthal correlation function from feed-down D-meson triggers, \({\tilde{C}}^{\mathrm{MC \,templ}}_{\mathrm{feed}-\mathrm{down}}(\Delta \varphi )\), was evaluated for each \(p_{\mathrm{T}} \) range with the PYTHIA6 event generator with Perugia-2011 tune at generator level (i.e. without detector effects and selection criteria). After being rescaled to the expected fraction of feed-down D-meson candidate triggers, \(1-f_{\mathrm{prompt}}\), this contribution was subtracted from the purity-corrected azimuthal-correlation function. The fully-corrected, per-trigger azimuthal correlation function of prompt D mesons and charged particles was obtained with this procedure, as summarised in Eq. 2

$$\begin{aligned}&\frac{1}{N_{\mathrm{D}}} \times \frac{\mathrm{d}N^{\mathrm{assoc}}}{\mathrm{d}\Delta \varphi }(\Delta \varphi ) = \frac{1}{f_{\mathrm{prompt}}} \nonumber \\&\quad \times \left[ b_{\mathrm{B}-\mathrm{bias}}(\Delta \varphi )\times p_{\mathrm{prim}}(\Delta \varphi )\times \tilde{C}_{\mathrm{{signal}}}(\Delta \varphi , \Delta \eta ) \right. \nonumber \\&\quad \left. - \left( 1-f_{\mathrm{prompt}}\right) \tilde{C}^{\mathrm{MC\ templ}}_{\mathrm{{feed}-\mathrm{down}}}(\Delta \varphi ) \right] . \end{aligned}$$
(2)

3.2 Details specific to the multiplicity-differential analysis

For the multiplicity-differential analysis, only \(\mathrm D^{0} \) mesons and their charge conjugates were used as trigger particles. The same selection criteria chosen for the multiplicity-integrated case were used for the \(\mathrm D^{0} \)-meson candidate triggers and the associated charged particles. The evaluation of the two-dimensional correlation function for the peak region and the sidebands of the \(\mathrm D^{0} \)-meson invariant-mass distributions was performed independently in each of the four V0M multiplicity classes, considering the same \(p_{\mathrm{T}} \) ranges of the multiplicity-integrated analysis, with the exception of \(24< p_{\mathrm{T}}^{\mathrm{D}} < 36\) GeV/\(c\), where the amount of collected data was not sufficient for performing the study. The invariant-mass distributions of the V0M multiplicity class I were characterised by a larger data sample and a larger statistical significance of the \(\mathrm D^{0} \) mass peak, profiting from the usage of the V0HM trigger, although they also showed a lower signal-to-background ratio due to the enhanced underlying-event activity.

For each V0M multiplicity class, the per-trigger azimuthal correlation function was obtained following Eqs. 1 and 2, but some of the quantities entering these equations were evaluated with a modified procedure. For the mixed-event correction, a different classification of the events in terms of SPD tracklet multiplicity needed to correlate \(\mathrm D^{0} \) mesons and charged particles from events with similar features was considered for each V0M multiplicity class, since the SPD tracklet multiplicity distributions obtained for each V0M multiplicity class were significantly different. The MC events used to evaluate the values of 1/(\(A\times \epsilon )^{\mathrm{trig}}\), for each V0M multiplicity class, were reweighted in order to reproduce the corresponding SPD tracklet multiplicity measured in data. The same values of charged-particle reconstruction and selection efficiency were instead used for the four V0M multiplicity classes since a negligible dependence of the efficiency on the event multiplicity was found in previous ALICE studies in the same collision system [44].

The purity correction, \(p_{\mathrm{prim}}(\Delta \varphi )\), was evaluated independently for each V0M multiplicity class, by applying the same MC reweighting procedure used for the \(\mathrm D^{0} \)-meson efficiency evaluation. A very small dependence on event multiplicity was obtained, with the overall differences smaller than \(0.5 \%\) between the values obtained in the four V0M multiplicity classes and those evaluated for the multiplicity-integrated analysis. Similarly, the correction factor \(b_{\mathrm{B}-\mathrm{bias}}(\Delta \varphi )\) was estimated separately for each of the four V0M multiplicity classes. A slight increase of the correction of about 1–2% depending on the \(p_{\mathrm{T}} \) range was obtained with decreasing event multiplicity. The largest value of \(b_{\mathrm{B}-\mathrm{bias}}(\Delta \varphi )\) was about \(5.5\%\) for V0M multiplicity class IV at \(\Delta \varphi =0\) for the lowest \(p_{\mathrm{T}}^{\mathrm{D}} \) interval and largest \(p_{\mathrm{T}}^{\mathrm{assoc}} \) interval. The feed-down subtraction procedure was left unaltered for the evaluation of the central values of the correlation function, assuming no modification of the prompt \(\mathrm D^{0} \)-meson fraction and beauty-quark fragmentation with event multiplicity. However, as described in more detail in Sect. 4, an additional systematic uncertainty was considered, accounting for a possible variation of the feed-down-to-prompt \(\mathrm D^{0} \)-meson production ratio with event multiplicity, which would impact the value of \(f_{\mathrm{prompt}}\).

3.3 Quantitative evaluation of correlation peak features

Based on their consistency within uncertainties, an average of the azimuthal-correlation functions for \(\mathrm D^{0} \), \(\mathrm D^{+} \), and \(\mathrm D^{*}(2010)^{+} \) meson triggers was evaluated for the multiplicity-integrated analysis. The average was obtained by weighting the correlation function from each species by the inverse of the quadratic sum of its statistical and systematic uncertainties uncorrelated among the three D-meson species. For each \(p_{\mathrm{T}} \) interval, the averaged azimuthal correlation function was fitted with the following function

$$\begin{aligned} f(\Delta \varphi )&= a + \frac{Y_{\mathrm{NS}}\times \beta }{2\alpha \Gamma (1/\beta )}\times e^{-\left( \frac{\Delta \varphi }{\alpha }\right) ^\beta } \nonumber \\&\qquad + \frac{Y_{\mathrm{AS}}}{\sqrt{2\pi }\sigma _{\mathrm{AS}}}\times e^{-\frac{(\Delta \varphi - \pi )^2}{2\sigma ^2_{\mathrm{AS}}}}. \end{aligned}$$
(3)

This function is composed of a generalised-Gaussian component for the description of the near-side peak (with the mean fixed at \(\Delta \varphi = 0\)), a Gaussian component for the away-side peak (with the mean fixed at \(\Delta \varphi = \pi \)), and a constant term (baseline) to account for the flat contribution that lies beneath the two correlation peaks. To grant sufficient stability to the fit, the \(\beta \) parameter of the generalised Gaussian was fixed to the value obtained for the correlation distribution predicted by PYTHIA8 simulations at generator level. In Eq. 3, the baseline value a was fixed to the weighted average of the points in the range \(\pi /4< |\Delta \varphi | < \pi /2\) (transverse region), to reduce the fit sensitivity to statistical fluctuations. The inverse of the squared statistical uncertainties of the points were used as weights. The fit to the correlation function allowed the extraction of quantitative observables that characterise the correlation peaks. In particular, the near- and away-side peak yields were obtained as the integral of the components describing each correlation peak, and their widths were parameterised by the quantities \(\alpha \sqrt{\Gamma (3/\beta )/\Gamma (1/\beta )}\) (square root of the generalised-Gaussian variance) and \(\sigma _{\mathrm{AS}}\), respectively.

4 Systematic uncertainties

4.1 Systematic uncertainties for the multiplicity-integrated analysis

The azimuthal correlation function obtained from the multiplicity-integrated analysis is affected by several systematic uncertainties due to the specific procedure and assumptions introduced for its evaluation. In the following, the approach used to estimate each systematic uncertainty source is briefly described.

The evaluation of \(S_{\mathrm {peak\ region}}\) and \(B_{\mathrm {peak\ region}}\) from the fit to the D-meson invariant-mass distributions (as described in Sect. 3.1 and Eq. 1) introduced a systematic uncertainty on the correlation function. The uncertainty was estimated by varying the fit procedure, specifically, by modelling the background distribution with a linear function or a second-order polynomial function instead of an exponential function (for \(\mathrm D^{0} \) and \(\mathrm D^{+} \) mesons only, where there is not a straightforward choice for the background fit function), considering a different histogram binning, varying the fit range, fixing the mean of the Gaussian term describing the mass peak to the world-average D-meson mass [41], or fixing the Gaussian width to the value obtained from MC studies. A systematic uncertainty ranging from 0.5 to 1.5%, depending on the \(p_{\mathrm{T}}^{\mathrm{D}} \) range and similar for all D-meson species, was estimated from the corresponding variation of the azimuthal-correlation function. No dependence on \(\Delta \varphi \) was found.

A 0.5–2% systematic uncertainty, depending on \(p_{\mathrm{T}}^{\mathrm{D}} \) range and D-meson species, was introduced due to the possible dependence of the shape of the background correlation function on the invariant-mass value of the trigger D meson. This source of uncertainty was estimated by evaluating \(\tilde{C}_{\mathrm{{sidebands}}}(\Delta \varphi , \Delta \eta )\) considering different invariant-mass sideband ranges. For \(\mathrm D^{0} \) and \(\mathrm D^{+} \) mesons, for which a sideband was defined on each side of the invariant-mass signal peak, \(\tilde{C}_\mathrm{{sidebands}}(\Delta \varphi , \Delta \eta )\) was also evaluated considering only the left or the right sideband. No azimuthal dependence was observed for this uncertainty.

The evaluation of the associated-particle reconstruction efficiency via MC studies introduced a further systematic uncertainty, estimated by varying the quality selection criteria applied on the reconstructed tracks, i.e. removing or tightening the request on minimum number of ITS clusters, requiring a hit on at least one of the two SPD layers, or varying the request on the minimum number of space points reconstructed in the TPC. An uncertainty of 4–5% was estimated, independent of the D-meson species, and no significant trend in \(\Delta \varphi \) was observed.

A systematic uncertainty affecting the D-meson reconstruction efficiency, related to potentially different distributions of the topological variables in MC and data, was estimated by testing a set of tighter and looser topological selections on D-meson candidates. An uncertainty ranging from 0.5 to 2% was assigned and the effect on the azimuthal correlation function was found to be \(\Delta \varphi \) independent.

The correlation function has an uncertainty related to the evaluation of the residual contamination from secondary particles [43]. To determine this, the analysis was repeated by varying the DCA selection in the xy plane from 0.1 cm to 2.4 cm and re-evaluating the purity correction of primary tracks for each variation. This resulted in a maximum, \(\Delta \varphi \)-independent, systematic uncertainty of 2% on the azimuthal-correlation function.

In addition to the above contributions, which all act as a scale uncertainty, the azimuthal correlation function is also affected by \(\Delta \varphi \)-dependent systematic uncertainties. The uncertainty on the evaluation of the beauty feed-down contribution to the azimuthal correlation function was determined by employing alternate templates of feed-down azimuthal-correlation functions, obtained from different event generators (PYTHIA6 with the Perugia-2010 tune [39] and PYTHIA8 with the 4C tune [45]), and by varying the value of \(f_\mathrm{prompt}\) following the procedure described in Ref. [42]. A \(\Delta \varphi \)-dependent uncertainty was obtained with a maximum value of 5%. The near-side region for the feed-down D-meson component of the correlation function was affected by a bias, favouring topologies with a small opening angle between the D meson and the other beauty-hadron decay products. This was corrected for as discussed in Sect. 3.1, and a \(\Delta \varphi \)-dependent bilateral and symmetric uncertainty for a possible over- or under-correction of this bias was evaluated as detailed in Ref. [21]. The largest value of the uncertainty was found to be 2.5% for \(\Delta \varphi \approx 0\).

The estimated systematic uncertainty values from each of the above sources affecting the azimuthal correlation function are summarised in Table 2. The overall systematic uncertainty in each \(\Delta \varphi \) bin of the correlation function was obtained as the sum in quadrature of the aforementioned contributions.

The systematic uncertainties on the peak observables were evaluated by considering several contributions: (i) the impact on the physical observables induced by the baseline position was estimated by considering alternate \(\Delta \varphi \) ranges for determining its value and repeating the fit; similarly, the possible bias induced by fixing the \(\beta \) parameter of the near-side to the predicted PYTHIA8 value was estimated by allowing \(\beta \) to vary within \(\pm 20\%\) of that value. For each observable, the root-mean-square of the relative variations from these alternative fits with respect to the central value of the observable was considered; (ii) the impact of the \(\Delta \varphi \)-dependent uncertainty on the correlation function was accounted for by coherently shifting its points to the upper and lower edges of their \(\Delta \varphi \)-dependent systematic uncertainty values. The fit was repeated and the variation of each observable value with respect to the default value was considered in each direction; (iii) the overall \(\Delta \varphi \)-independent systematic uncertainty acts as a scaling factor on the correlation function, hence it impacted the near- and away-side peak yield values by the same relative amount. The overall systematic uncertainty on the peak yields was therefore obtained by summing in quadrature the contributions from (i), (ii), and (iii). For the near-side and away-side widths, which are insensitive to scale factors, the sum in quadrature of only the contributions from (i) and (ii) was considered.

4.2 Systematic uncertainties for the multiplicity-differential analysis

Some of the systematic uncertainties affecting the correlation function were estimated separately in each of the V0M multiplicity class, following the same prescriptions described in Sect. 4.1. In particular, this was done for the uncertainties on the yield extraction and background correlation shape since they are related to the features of the invariant-mass distributions, which show a significant multiplicity evolution. The uncertainty related to the bias affecting the topological selection of feed-down \(\mathrm D^{0} \) mesons was also re-evaluated, since a slight multiplicity dependence of the related correction was found. Similar values of the uncertainty were obtained compared to the multiplicity-integrated case. For the subtraction of the beauty feed-down contribution, an additional systematic uncertainty was considered, related to a possible multiplicity dependence of the relative fraction of feed-down \(\mathrm D^{0} \) mesons in the \(\mathrm D^{0} \)-meson raw yields that determines the amount of the feed-down contribution. The evaluation of this uncertainty followed a similar procedure as the one described in Ref. [9], and led to an asymmetry of the feed-down systematic uncertainty, which increased up to \(^{+5\%}_{-9\%}\) for the V0M multiplicity class I. For the other systematic uncertainty sources affecting the azimuthal-correlation function, the same values estimated for the multiplicity-integrated analysis were adopted. The uncertainty values for the multiplicity-differential analysis are reported in Table 2.

Table 2 The list of the systematic uncertainty contributions affecting the azimuthal correlation function and their typical values. If not specified, the uncertainties do not depend on \(\Delta \varphi \)
Full size table

The evaluation of the systematic uncertainties on the near- and away-side peak observables was unmodified with respect to the procedure described in Sect. 4.1 for the multiplicity-integrated analysis. In addition, the impact on the peak observables related to the possible presence of long-range azimuthal correlations between \(\mathrm D^{0} \) mesons and charged particles in high-multiplicity collisions was studied by replacing in Eq. 3 the constant term with a \(v_{2\Delta }\)-like modulation. The elliptic-flow coefficient values of \(\mathrm D^{0} \) mesons and charged particles adopted to evaluate the modulation were defined based on the measurements in Ref. [10] and Ref. [46], respectively. For the V0M multiplicity class I, variations within 2% were found for all observables and \(p_{\mathrm{T}} \) ranges except for the range 3 \(< p_{\mathrm{T}}^{\mathrm{D}}<\) 5 GeV/\(c\), where a reduction of the near- and away-side peak yields as large as 8% was observed. These differences were assigned as a systematic uncertainty.

5 Results

5.1 Multiplicity-integrated results in pp collisions at \(\sqrt{s} = 13\) TeV

For the multiplicity-integrated analysis, the azimuthal correlation function of D mesons with charged particles was computed for the five D-meson \(p_{\mathrm{T}} \) intervals \(3< p_{\mathrm{T}}^{\mathrm{D}} < 5\) GeV/\(c\), \(5< p_{\mathrm{T}}^{\mathrm{D}} < 8\) GeV/\(c\), \(8<p_{\mathrm{T}}^{\mathrm{D}} < 16\) GeV/\(c\), \(16< p_{\mathrm{T}}^{\mathrm{D}} < 24\) GeV/\(c\), and \(24< p_{\mathrm{T}}^{\mathrm{D}} < 36\) GeV/\(c\), and for associated particle \(p_{\mathrm{T}} \) range \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\) and the sub-intervals \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), \(1<p_{\mathrm{T}}^{\mathrm{assoc}} < 2\) GeV/\(c\), and \(2<p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\).

Fig. 1
figure 1

Average of the azimuthal-correlation functions of \(\mathrm D^{0} \), \(\mathrm D^{+} \), and \(\mathrm D^{*+} \) mesons with associated particles, after the baseline subtraction, in pp collisions at \(\sqrt{s}\) = 5.02 [21], 7 [3], and 13 TeV, for 3 \(< p_{\mathrm{T}}^{\mathrm{D}}<\) 5 GeV/\(c\), 8 \(< p_{\mathrm{T}}^{\mathrm{D}}<\) 16 GeV/\(c\), and 16 \(<p_{\mathrm{T}}^{\mathrm{D}}<\) 24 GeV/\(c\) (from left to right) and 0.3 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 1 GeV/\(c\), 1 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 2 GeV/\(c\) (top and bottom panels, respectively). Data at \(\sqrt{s}\) = 7 TeV are not available for all the \(p_{\mathrm {T}}\) regions. Statistical and \(\Delta \varphi \)-dependent systematic uncertainties are shown as vertical error bars and boxes, respectively, and \(\Delta \varphi \)-independent uncertainties are written as text. The uncertainties from the subtraction of the baseline are displayed as boxes at \(\Delta \varphi >\pi \)

Full size image
Fig. 2
figure 2

Near-side (left panel) and away-side (right panel) peak yields (first row) and widths (third row) obtained from a fit to the azimuthal correlation function after the baseline subtraction. The measurements are compared with ALICE results obtained in pp collisions at \(\sqrt{s}\) = 5.02 TeV [21] and 7 TeV [3], for 0.3 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 1 GeV/c, 1 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 2 GeV/c. Only near-side observables were computed in the \(\sqrt{s}\) = 7 TeV analysis. The ratios to the \(\sqrt{s} = 5.02\) TeV results are shown in the second and fourth rows for yields and widths, respectively

Full size image

Figure 1 shows examples of correlation functions obtained from the analysis in pp collisions at a centre-of-mass energy of \(\sqrt{s} \) = 13 TeV, compared with results previously reported by ALICE in pp collisions at \(\sqrt{s} \) = 5.02 TeV [21] and \(\sqrt{s} \) = 7 TeV [3] (the latter is available only for two kinematic ranges). The baseline value is closely related to the number of charged particles produced at midrapidity, and therefore has a strong dependence on \(\sqrt{s} \), due to the increase of charged-particle production with increasing centre-of-mass energy [47]. It was subtracted from the correlation functions in order to focus the comparison on the peak features. The shape of the distribution after the baseline subtraction and the properties of the correlation peaks at the three centre-of-mass energies agree within uncertainties. The analysis at \(\sqrt{s} = 13\) TeV profits from a larger data sample, resulting in substantially smaller point-by-point statistical fluctuations (up to 50% with respect to \(\sqrt{s} = 5.02\) TeV results in the range \(16<p_{\mathrm{T}}^{\mathrm{D}} <24\) GeV/\(c\)), and leading to smaller uncertainties from the subtraction of the baseline.

More quantitative results are provided by the comparison of the near- and away-side peak yields and widths in pp collisions at different centre-of-mass energies, presented in Fig. 2. The results for pp collisions at \(\sqrt{s} = 7\) TeV were obtained after refitting the correlation functions measured in Ref. [3] with the improved function described in Eq. 3, and evaluating the systematic uncertainties accordingly. The near-side peak yield values obtained from the \(\sqrt{s} = 13\) TeV data are compatible within the uncertainties with those at lower energies, exhibiting the same increasing trend of the yields with \(p_{\mathrm{T}}^{\mathrm{D}} \) in both the \(p_{\mathrm{T}}^{\mathrm{assoc}} \) intervals shown. An overall agreement is also observed between the \(\sqrt{s} = 5.02\), 7, and 13 TeV results for the near-side widths. From the \(\sqrt{s} = 13\) TeV results, an indication of a near-side peak narrowing for increasing \(p_{\mathrm{T}}^{\mathrm{D}} \) emerges, in particular for the \(1< p_{\mathrm{T}}^{\mathrm{assoc}} < 2\) GeV/\(c\) range, that was not observed at lower energies because of the lower precision. This peak narrowing could be originated by two simultaneous effects: (i) a more collimated angular pattern of the partons fragmented from charm quarks; (ii) an increased collinearity of charm and anti-charm quarks produced from gluon-splitting mechanism. Both effects are related to the increased boost, on average, of the fragmenting (splitting) parton when considering D-meson triggers with larger \(p_{\mathrm{T}} \). An agreement within the uncertainties is also observed for the away-side peak results, and similar conclusions as those expressed for the near-side peak can be drawn. For the away-side observables only results at \(\sqrt{s} = 5.02\) TeV are available for the comparison (and only for a restricted set of kinematic ranges), since the azimuthal-correlation functions for pp collisions at \(\sqrt{s} = 7\) TeV were not precise enough to allow the characterisation of the away-side region, as discussed in Ref. [3].

A similar comparison was performed using simulated pp collisions obtained with PYTHIA6 (with Perugia-2011 tune) and POWHEG+PYTHIA6 [19, 20, 48] event generators. A slight increase of the near-side yield values (5–10% depending on the \(p_{\mathrm{T}} \) range) and a mild decrease of the away-side yield values (10–15%) was observed when increasing the centre-of-mass energy from \(\sqrt{s} = 5.02\) to 13 TeV, with small differences between the two generators. This could be ascribed to an increased contribution of NLO production processes (already included in the hard scattering in POWHEG+PYTHIA6, and accounted for during the parton-shower development in PYTHIA6), as well as to a harder charm-quark \(p_{\mathrm{T}} \) spectrum at larger centre-of-mass energies. These differences are within the overall precision of the data measurements. No visible energy dependence for both near- and away-side peak widths was found, for both generators, as observed in data.

5.2 Results for different V0M multiplicity classes

The azimuthal-correlation functions evaluated in the four classes of V0M multiplicity are compared in Fig. 3, for the four available \(\mathrm D^{0} \)-meson \(p_{\mathrm{T}} \) ranges (one per column) and the four different \(p_{\mathrm{T}}^{\mathrm{assoc}} \) intervals (one per row). The baseline value largely increased from V0M multiplicity class IV towards V0M multiplicity class I, as expected due to the very different underlying-event activity, and it was subtracted from the correlation functions. This comparison suggests a similar shape and \(p_{\mathrm{T}} \) evolution of the azimuthal-correlation function, as well as consistent near- and away-side features for the four V0M multiplicity classes.

Fig. 3
figure 3

Azimuthal-correlation functions of \(\mathrm D^{0} \) mesons with associated particles, after the subtraction of the baseline, in pp collisions at \(\sqrt{s}\) = 13 TeV, for 3 \(< p_{\mathrm{T}}^{\mathrm{D}}<\) 5 GeV/\(c\), 5 \(< p_{\mathrm{T}}^{\mathrm{D}}<\) 8 GeV/\(c\), 8 \(< p_{\mathrm{T}}^{\mathrm{D}}<\) 16 GeV/\(c\), and 16 \(<p_{\mathrm{T}}^{\mathrm{D}}<\) 24 GeV/\(c\) (from left to right) and \(p_{\mathrm{T}}^{\mathrm{assoc}}>\) 0.3 GeV/\(c\), 0.3 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 1 GeV/\(c\), 1 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 2 GeV/\(c\), 2 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 3 GeV/\(c\) (from top to bottom) for different multiplicity classes estimated with V0M. The four multiplicity classes are shown with different marker styles. Statistical and \(\Delta \varphi \)-dependent systematic uncertainties are shown as vertical error bars and boxes, respectively, and \(\Delta \varphi \)-independent uncertainties are written as text. The uncertainties from the subtraction of the baseline are displayed as boxes at \(\Delta \varphi >\pi \)

Full size image

The near-side peak yield and width values obtained from the fit to the azimuthal correlation function in the four V0M multiplicity classes, for the same kinematic ranges, are shown in Fig. 4 (first and third rows, respectively). Apart from a tension for low \(p_{\mathrm{T}}^{\mathrm{D}} \), 2 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 3 GeV/\(c\), the yield measurements follow a similar increasing trend with \(p_{\mathrm{T}}^{\mathrm{D}} \). Similar values are observed from the near-side peak widths. These results indicate no significant modification of the charm fragmentation and hadronisation in collisions of varying charged-particle multiplicities. The ratios of the yield and width results in V0M classes II, III, and IV over those in V0M class I, shown in the second and fourth rows of Fig. 4, respectively, also confirm this conclusion.

The evaluation of away-side peak yields and widths as a function of the event multiplicity were performed only in the integrated associated particle interval \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\), due to their large sensitivity to point-by-point statistical fluctuations. The away-side peak observable values for each of the four V0M multiplicity classes are shown in Fig. 5, together with the ratios of the values in the V0M classes II, III and IV over the values in the V0M class I. As observed for the near-side, the same increasing trend with \(\mathrm D^{0} \) \(p_{\mathrm{T}} \) is present for the yields, among the four V0M multiplicity classes. A hint of narrowing of the away-side peak, visible in the multiplicity-integrated results at \(\sqrt{s} = 13\) TeV, can also be seen in the multiplicity-dependent results. The away-side yield and width values are fully consistent within the uncertainties among all four V0M classes.

Though with sizeable uncertainties, these measurements point towards consistency of the jet-induced correlation peak structure and shape in high- and low-multiplicity events, and thus contribute to confirm the assumptions done in the measurements of the elliptic-flow coefficient of charm particles in high-multiplicity pp collisions [10, 11].

Fig. 4
figure 4

Near-side associated peak yields (top row) and widths (third row) measured in pp collisions at \(\sqrt{s}\) = 13 TeV, for the four V0M multiplicity classes, shown with different marker styles. The ratios of yield (width) values in each V0M class with respect to those in the V0M class I are shown in the second (fourth) row. Results are presented as a function of the \(\mathrm D^{0} \)-meson \(p_{\mathrm{T}} \), for \(p_{\mathrm{T}}^{\mathrm{assoc}}>\) 0.3 GeV/\(c\) and the sub-ranges 0.3 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 1 GeV/\(c\), 1 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 2 GeV/\(c\), and 2 \(<p_{\mathrm{T}}^{\mathrm{assoc}}<\) 3 GeV/\(c\) (from left to right). Statistical and systematic uncertainties are shown as vertical error bars and boxes, respectively

Full size image
Fig. 5
figure 5

Away-side associated peak yields (left) and widths (right) measured in pp collisions at \(\sqrt{s}\) = 13 TeV, for the four V0M multiplicity classes, shown with different marker styles. The ratios of the observable values in each V0M class with respect to those in V0M class I are represented in the bottom insets. Results are presented as a function of the \(\mathrm D^{0} \)-meson \(p_{\mathrm{T}} \), for \(p_{\mathrm{T}}^{\mathrm{assoc}}>\) 0.3 GeV/\(c\). Statistical and systematic uncertainties are shown as vertical error bars and boxes, respectively

Full size image

5.3 Comparison of the ALICE results with model predictions

The near- and away-side peak yields measured in pp collisions at \(\sqrt{s} = 13\) TeV and reported in Sect. 5.1 were compared to predictions from several event generators. This allowed for verifying, for each model, whether its specific implementation of the processes leading from charm-quark production to final-state particles was adequate for describing the measured observables. A detailed description of the models used for the comparison is provided in Ref. [21]. These models include PYTHIA8 [13] with 4C [45] tune, PYTHIA6 [12] with Perugia-2011 tune [39], POWHEG+PYTHIA8 [19, 20, 48] with hard-scattering matrix elements evaluated at NLO or at LO accuracy, HERWIG 7 [16, 17], and EPOS 3.117 [14, 15].

For each model, the average of the \(\mathrm D^{0} \), \(\mathrm D^{+} \), and \(\mathrm D^{*+} \) azimuthal-correlation functions with charged particles was evaluated, using the same prescriptions applied for data analysis in terms of kinematic and particle-species selections. The evaluation of the peak observables from the fit to the correlation distribution followed the same approach employed on data, except for the estimation of the baseline. Since the statistical fluctuations in the transverse region are negligible for the models, the minimum of the azimuthal correlation function was directly considered as the baseline value. A systematic uncertainty on the peak observables was then assigned by performing an alternate fit, fixing the baseline as the weighted average of the two lowest points of the azimuthal-correlation function.

The near-side peak observable trends for both models and data are illustrated in Fig. 6 as a function of the trigger D-meson \(p_{\mathrm{T}} \), in the \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\) interval and in the other three kinematic sub-ranges. The first and third rows show the yield and the width values, while the second and the fourth show the ratios of model predictions with respect to data. In these ratio panels, model statistical and systematic uncertainties are shown as error bars and boxes, respectively, while data statistical and systematic uncertainties are summed in quadrature, and the resulting uncertainty is represented as a solid grey band. The increasing trend of the near-side yield with the trigger particle \(p_{\mathrm{T}} \) seen in data is obtained by all the MC predictions, but with different strengths. A hierarchy can be observed for the yield values, with EPOS systematically providing the largest yields, followed by POWHEG+PYTHIA8 NLO, POWHEG+PYTHIA8 LO, and then by PYTHIA6 and PYTHIA8. HERWIG predicts the lowest yields for \(p_{\mathrm{T}}^{\mathrm{D}} < 8\) GeV/\(c\) and \(p_{\mathrm{T}}^{\mathrm{assoc}} > 1\) GeV/\(c\). Its predicted yield values for the other kinematic ranges, instead, are generally in between PYTHIA and POWHEG+PYTHIA8. The best description of the measurements is provided by the POWHEG+PYTHIA8 and by PYTHIA generators, with POWHEG+PYTHIA8 (both NLO and LO) generally performing better at lower \(p_{\mathrm{T}}^{\mathrm{D}} \) and PYTHIA (both versions) at higher \(p_{\mathrm{T}}^{\mathrm{D}} \). A slight difference is observed between the NLO and LO implementations of POWHEG+PYTHIA8, with the former providing larger yields (by 5 to 15%, increasing with the D-meson \(p_{\mathrm{T}})\), and overall providing a better description of the data in the lower \(p_{\mathrm{T}}^{\mathrm{D}} \) region, while the latter has a better agreement with data above 8 GeV/\(c\). These differences can be understood in terms of the different relative contribution of the NLO production mechanisms, as already discussed in Ref. [21]. HERWIG predictions tend to underestimate the value of the near-side yield in the kinematic region \(p_{\mathrm{T}}^{\mathrm{D}} < 8\) GeV/\(c\) and \(p_{\mathrm{T}}^{\mathrm{assoc}} > 1\) GeV/\(c\), while for the other kinematic regions the predictions are compatible with the data. EPOS predictions (not available for the range \(24< p_{\mathrm{T}} < 36\) GeV/\(c\)) overestimate the near-side yield measurements by a factor of about 2 through all the studied kinematic ranges. Generally, smaller differences are obtained for the widths of the various models with respect to those observed for the yields. POWHEG+PYTHIA8 NLO predicts the broadest near-side peaks. The near-side width data measurements hint towards a slight sharpening of the near-side peak width with increasing \(p_{\mathrm{T}}^{\mathrm{D}} \), while most of the models describe the width as nearly flat. However, all models are able to reproduce the measured width within the uncertainties.

Fig. 6
figure 6

Near-side associated peak yields (top row) and widths (third row down) in pp collisions at \(\sqrt{s}\) = 13 TeV, compared with predictions by PYTHIA, POWHEG+PYTHIA8, HERWIG, and EPOS 3 event generators with various configurations. Results are presented as a function of the D-meson \(p_{\mathrm{T}} \), for \(p_{\mathrm{T}}^{\mathrm{assoc}}>\) 0.3 GeV/\(c\), \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), \(1< p_{\mathrm{T}}^{\mathrm{assoc}} < 2\) GeV/\(c\), and \(2< p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\) (from left to right). The ratios of model predictions to data measurements for yield (width) values are shown in the second (fourth) row down. In these rows, model statistical and systematic uncertainties are shown as vertical error bars and boxes, respectively, while data total uncertainties are displayed as a solid band

Full size image
Fig. 7
figure 7

Away-side associated peak yields (top row) and widths (third row down) in pp collisions at \(\sqrt{s}\) = 13 TeV, compared with predictions by PYTHIA, POWHEG+PYTHIA8, HERWIG, and EPOS 3 event generators with various configurations. Results are presented as a function of the D-meson \(p_{\mathrm{T}} \), for \(p_{\mathrm{T}}^{\mathrm{assoc}}>\) 0.3 GeV/\(c\), \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), \(1< p_{\mathrm{T}}^{\mathrm{assoc}} < 2\) GeV/\(c\), and \(2< p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\) (from left to right). The ratios of model predictions to data measurements for yield (width) values are shown in the second (fourth) row down. In these rows, model statistical and systematic uncertainties are shown as vertical error bars and boxes, respectively, while data total uncertainties are displayed as a solid band

Full size image

A similar comparison for the away-side peak is shown in Fig. 7. POWHEG+PYTHIA8 NLO and LO implementations provide the highest away-side yields, with the LO implementation generally 5% above the NLO one, possibly due to an increased amount of back-to-back production processes. PYTHIA6 and PYTHIA8 generally provide slightly lower away-side yield values than POWHEG+PYTHIA8 NLO and LO expectations, with PYTHIA8 tending to be on the lower side compared to PYTHIA6. HERWIG predicts values lower than all the other models (about 20% lower than POWHEG+PYTHIA8 NLO yields). EPOS predicts a stronger increasing trend of the away-side yields with \(p_{\mathrm{T}}^{\mathrm{D}} \) than the other generators. For this observable, data uncertainties are not small enough to be sensitive to all the differences highlighted above. The best agreement with data is provided by PYTHIA6, PYTHIA8, and HERWIG. The POWHEG+PYTHIA8 predicted yields, for both implementations, are on the higher side compared to data, in particular for lower associated-particle \(p_{\mathrm{T}} \). EPOS predictions tend to underestimate the yield for \(p_{\mathrm{T}}^{\mathrm{D}} < 5\) GeV/\(c\), while for \(16< p_{\mathrm{T}}^{\mathrm{D}} < 24\) GeV/\(c\) its away-side yield predictions are higher than the measured values. The narrowing of the away-side peak with increasing \(p_{\mathrm{T}}^{\mathrm{D}} \) is clearly evident in the third row of Fig. 7, both for model predictions and data, except for \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), where all models predict a rather flat trend. Only EPOS has a slightly different behaviour compared to the other models, showing an approximately flat \(p_{\mathrm{T}}^{\mathrm{D}} \) trend of the away-side width over the full \(p_{\mathrm{T}}^{\mathrm{assoc}} \) range, albeit with quite large model uncertainties. In terms of absolute values, PYTHIA6 away-side width expectations are systematically higher than all the other models, by an overall 20\(\%\), with increasing differences for increasing \(p_{\mathrm{T}}^{\mathrm{assoc}} \), and also tend to overestimate the measured width values. All the other models predict similar values, consistent with data except for \(p_{\mathrm{T}}^{\mathrm{D}} > 16\) GeV/\(c\), where they slightly overestimate the data measurements.

Fig. 8
figure 8

Near- and away-side peak yields (first column), widths (second column), near-side peak \(\beta \) parameter and baseline (third column) from fits to the D-meson and charged particle azimuthal-correlation function, from PYTHIA8 simulations obtained with different parton-level contributions. The predictions are obtained for multiplicity-integrated pp collisions at \(\sqrt{s} = 13\) TeV, as a function of the D-meson \(p_{\mathrm{T}} \), for \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\), and compared with ALICE data

Full size image
Fig. 9
figure 9

Near- and away-side peak yields (first column), widths (second column), near-side peak \(\beta \) parameter and baseline (third column) from fits to the D-meson and charged particle azimuthal-correlation function, from POWHEG+PYTHIA8 simulations obtained with different parton-level contributions. The predictions are obtained for multiplicity-integrated in pp collisions at \(\sqrt{s} = 13\) TeV, as a function of the D-meson \(p_{\mathrm{T}} \), for \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\), and compared with ALICE data

Full size image

6 Parton-level studies with PYTHIA8 and POWHEG+PYTHIA8

From the comparative studies discussed in Sect. 5.3, PYTHIA8 and POWHEG+PYTHIA8 provide, overall, the most accurate description of the near- and away-side correlation peak features. A more detailed investigation was performed using these models to expose how different stages of the event generation that determine the formation of the final-state particles generally influence the development of the features of the correlation peak and the azimuthal-correlation function.

At large momentum or short distances, e.g. the hard-parton scattering leading to heavy quark production, QCD is asymptotically free. It implies that the coupling constant is small, so a perturbative approach can be employed. Before hard scattering takes place, partons from the incident beam protons can radiate gluons in the so-called initial-state radiation (ISR) process. Similarly, outgoing partons from the hard scatterings can produce a shower of softer particles via a final-state radiation (FSR) process. Since hadrons are composite objects, more than one distinct hard-parton interaction can occur in a pp collision, and proton remnants can also scatter again on each other. Such processes are called multi-parton interactions (MPI), and are responsible for the production of a large fraction of the particles uncorrelated with the D-meson candidate trigger, giving a substantial contribution to the underlying event (UE) observed in the collision final state. Additionally, as detailed in Ref. [9], in the MPI implementation used in PYTHIA8 (which also drives the MPI process in POWHEG+PYTHIA8 simulations) charm-quark production can occur not only from the first (hardest) hard scattering, but also from hard processes in the various MPI occurring in the collisions, ordered with decreasing hardness. There is also some correlation between FSR+ISR and MPI processes, since initial- and final-state radiations are generated from all the parton interactions occurring in the collision, and are thus enhanced in presence of MPI, further increasing the collision multiplicity.

In the following, PYTHIA8 and POWHEG+PYTHIA8 predictions, for standard simulations and for events generated after deactivating the previously mentioned parton-level contributions (FSR+ISR, and FSR+ISR+MPI), are compared. The latter case, in particular, provides a detailed view of the correlation function from the hard-scattering outgoing partons, though hadronisation and decays are still present and partially modify the original, parton-level, correlation shape. The same procedure was performed for the predictions of near- and away-side peak yields, widths, and baseline value, which are then compared with the data, and of the \(\beta \) parameter of the near-side peak.

Figure 8 shows the observables extracted from the fit to the azimuthal correlation function for different parton-level contributions from PYTHIA8 event generator, compared with the data, for \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\) and as a function of \(p_{\mathrm{T}}^{\mathrm{D}} \). In the first column, at high \(p_{\mathrm{T}}^{\mathrm{D}} \) the near-side and away-side yields show no relevant contribution of MPI, while FSR and ISR processes contribute to an increase of both peak yields, as expected, since for high-momentum partons they lead to additional production of collinear particles. Even at high \(p_{\mathrm{T}}^{\mathrm{D}} \), with all three processes switched off, PYTHIA8 predicts peak yields amounting to roughly half of the measured yields. In the lowest \(p_{\mathrm{T}}^{\mathrm{D}} \) interval, instead, FSR, ISR, and MPI lead to a decrease of the peak yields. In the second column, the near-side width shows no significant modification for the various configurations, apart from a slight increase of the width observed when switching off MPI. These insights point towards a relevant role of hadronisation, which remains in place also in the absence of FSR, ISR, and MPI processes, in shaping the near-side correlation peak. As expected, the away-side peak is wider than the near-side peak because of a combined contribution of parton-level effects. When FSR and ISR are turned off, the peak width is substantially decreased. This could be explained by the lack of radiation (in particular hard gluon emissions) decreasing the deflection angle of the recoil jet. In the third column, a mild dependence is observed for the near-side \(\beta \) parameter. In the simulations, the available sample is much larger than in data, so the \(\beta \) parameter could be left free in the fit function, and its value is compared among the different model configurations. A very small contribution of parton-level processes to the baseline comes directly from the hard scattering (and subsequent hadronisation and decays), as expected since PYTHIA8 is a LO generator. All the other processes (ISR, FSR, and MPI) contribute to the baseline, with the MPI roughly equivalent to the sum of ISR and FSR. This is expected since MPI mostly affect the underlying event, whose particles point to directions largely unrelated to the trigger D-meson one. The further increase of the baseline when activating ISR and FSR processes is due to the fact that, as mentioned above, they act on all parton scatterings, including those occurring in MPI.

Figure 9 shows the same set of observables as shown in Fig. 8 and provides the same comparison to data, but with predictions obtained from POWHEG+PYTHIA8. Compared to Fig. 8, a larger near-side peak yield is already obtained when ISR, FSR, and MPI processes are switched off. This can be explained with the inclusion of NLO processes directly in the hard scattering in POWHEG, rather than reproducing their effect in the parton shower as in PYTHIA8. Also in this case, MPI does not contribute to the peak yield. A similar behaviour as that of Fig. 8 is observed for the near-side peak width. The near-side peak \(\beta \) parameter value shows a decreasing influence of ISR, FSR, and MPI processes for increasing \(p_{\mathrm{T}}^{\mathrm{D}} \). For the away-side peak, no major differences with respect to PYTHIA8 are found for the yield contributions, while a slightly lower influence of ISR and FSR is obtained for the widths, which can also be understood with the above consideration. For the baseline, a higher value with respect to PYTHIA8 is obtained when switching off all the parton-level processes, consistent with the non back-to-back topology of NLO processes (in particular, flavour excitation with a nearly flat contribution in \(\Delta \varphi \) [1]). The FSR and ISR contributions to the baseline increase with \(p_{\mathrm{T}}^{\mathrm{D}} \), in contrast to PYTHIA8, and partially drive the rising \(p_{\mathrm{T}}^{\mathrm{D}} \) trend of the baseline observed for the full simulation.

The comparison of the contributions of the various processes to the correlation function helps in understanding better the source of the difference between PYTHIA8 and POWHEG+PYTHIA8 simulations observed in Sect. 5.3. Figure 10 displays the azimuthal correlation function of D mesons with charged particles obtained from PYTHIA8 and POWHEG+PYTHIA8 simulations sequentially deactivating the different parton-level contributions in pp collisions at \(\sqrt{s} = 13\) TeV for \(0.3<p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\) and \(3<p_{\mathrm{T}}^{\mathrm{D}} < 5\) GeV/\(c\). Figure 11 shows the same quantities for a higher momentum range, i.e. \(2<p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\) and \(24<p_{\mathrm{T}}^{\mathrm{D}} < 36\) GeV/\(c\).

Fig. 10
figure 10

Azimuthal correlation function of D mesons with charged particles with different parton-level contributions from PYTHIA8 (left panel) and POWHEG+PYTHIA8 simulations (right panel) in pp collisions at \(\sqrt{s} = 13\) TeV, for \(0.3<p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\) and \(3<p_{\mathrm{T}}^{\mathrm{D}} < 5\) GeV/\(c\)

Full size image
Fig. 11
figure 11

Azimuthal correlation function of D mesons with charged particles with different parton-level contributions from PYTHIA8 (left panel) and POWHEG+PYTHIA8 simulations (right panel) in pp collisions at \(\sqrt{s} = 13\) TeV, for \(2<p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\) and \(24<p_{\mathrm{T}}^{\mathrm{D}} < 36\) GeV/\(c\)

Full size image

Most of the features expressed above are seen in a more qualitative but clearer way in these figures. In particular, the larger baseline in the case without FSR, ISR, and MPI for POWHEG+PYTHIA8 is clearly visible for the low-\(p_{\mathrm{T}} \) range. In general, for both simulations the relevant ISR, FSR, and MPI contributions to the baseline are obtained over the whole \(\Delta \varphi \) range when focusing in the low-momentum region, while most of the off-peak contribution disappears when considering higher-\(p_{\mathrm{T}} \) trigger and associated particles. In the high-\(p_{\mathrm{T}} \) region, the difference in the contributions to the peaks between the generators is evident: in particular, for POWHEG+PYTHIA8 the near-side peak yield is already larger for the case without ISR, FSR, and MPI processes. Also, at high \(p_{\mathrm{T}} \) two sharp peaks appear for PYTHIA8, when parton showering and MPI effects are turned off, while for POWHEG+PYTHIA8 wider peaks emerge already from the hard-scattering, due to higher-order charm-production processes. The addition of the parton showering processes in PYTHIA8 allows to reconcile most of the differences in the correlation peak shape, in particular for the widths, while some residual differences remain present for the yields for the full simulation, as discussed in Sect. 5.3.

7 Summary

A study of the azimuthal correlation function of D mesons with charged particles, measured in pp collisions at \(\sqrt{s} = 13\) TeV with the ALICE detector, was presented. The pattern of the correlation function and the features of near- and away-side correlation peaks, extracted via a fit to the correlation function, were characterised in five D-meson momentum ranges, from 3 to 36 GeV/\(c\), and for the associated-particle range \(p_{\mathrm{T}}^{\mathrm{assoc}} > 0.3\) GeV/\(c\) and the three sub-ranges \(0.3< p_{\mathrm{T}}^{\mathrm{assoc}} < 1\) GeV/\(c\), \(1<p_{\mathrm{T}}^{\mathrm{assoc}} < 2\) GeV/\(c\), and \(2<p_{\mathrm{T}}^{\mathrm{assoc}} < 3\) GeV/\(c\).

The measurement precision is significantly improved compared to previous ALICE results in pp collisions at \(\sqrt{s} = 7\) TeV [3] and \(\sqrt{s} = 5.02\) TeV [21]. The correlation function shape, as well as the peak yields and widths, are compatible within uncertainties with those lower-energy measurements, confirming the expectation from PYTHIA8 and POWHEG+PYTHIA8 generators of little dependence on the collision energies.

The possible evolution of the correlation function with the event multiplicity was probed by performing the analysis in four multiplicity ranges, measured with the V0M estimator, profiting from a dedicated high-multiplicity trigger provided by the V0 detector. Though the uncertainties do not allow for firm conclusions, a strong variation of the correlation function with multiplicity is excluded, suggesting that when the charm quarks hadronise into \(\mathrm D^{0} \) mesons, the charm-quark fragmentation and hadronisation processes are not particularly sensitive to the event multiplicity. The overall compatibility of the correlation-peak features for different event multiplicities tends to support the scenario that, independently of the number of charm quarks produced in the collision, they undergo similar fragmentation and hadronisation into \(\mathrm D^{0} \) mesons. Such a scenario indicates an increased, but independent, charm production from MPI in high-multiplicity pp collisions, and is one of the mechanisms proposed to explain the observed trend of D-meson self-normalised yields at different relative multiplicities [9].

The measured near- and away-side peak yields and widths were compared, in the accessible kinematic ranges, with expectations from state-of-the-art models capable of producing peak observables, such as PYTHIA, POWHEG+PYTHIA, HERWIG, and EPOS. Among these models, PYTHIA8 and POWHEG+PYTHIA8 provide the best overall description of the data, HERWIG underestimates the near-side yields at low \(p_{\mathrm{T}}^{\mathrm{D}} \) and at high \(p_{\mathrm{T}}^{\mathrm{assoc}} \), while EPOS overestimates the near-side yields over the whole kinematic range, predicting also a steeper \(p_{\mathrm{T}}^{\mathrm{D}} \) dependence of the away-side yields. A study of the influence of parton-level processes, as ISR, FSR, and MPI, on the shape of the final-state correlation distribution and on its peak properties was also performed with PYTHIA8 and POWHEG+PYTHIA8. These studies are important not only for a better understanding of the underlying physics in pp collisions, but also for interpreting possible modifications of the correlation peaks in Pb–Pb collisions due to interactions of the heavy quarks in the quark–gluon plasma. This measurement is expected to become accessible during the LHC Run 3, with the upgraded ALICE detector [49, 50].