1 Introduction

In high-energy hadronic collisions, heavy quarks (charm and beauty) are mainly produced in hard parton scattering processes. Due to the large momentum transfer characterizing these processes, their inclusive production cross sections can be calculated in the framework of perturbative quantum chromodynamics (pQCD) [1,2,3,4,5]. The production cross sections of several open heavy-flavor hadrons and of their decay leptons in pp collisions were measured at both mid- and forward-rapidity at the LHC [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27], and are described by pQCD calculations [28,29,30] with large theoretical uncertainties. The charm-hadron production cross section calculations in the pQCD frameworks are based on the factorization of parton distribution functions (PDF), the partonic cross section, and the fragmentation function. Recent measurements of charm-baryon production at midrapidity in pp collisions [31,32,33,34,35,36,37,38,39,40,41,42] are not reproduced by pQCD calculations and event generators adopting a fragmentation model tuned on \(\textrm{e}^{+}\textrm{e}^{-}\) data. A better description of these measurements can be obtained by models including hadronization mechanisms such as quark coalescence [43], additional color reconnections among parton fragments [44], or by including enhanced feed-down from higher-mass charm-baryon states within a statistical hadronization approach [45], where the higher-mass excited charm-baryon states are predicted by the Relativistic Quark Model [46] but not yet measured. More differential measurements are needed to better understand the fragmentation (parton showering) and hadronization of heavy quarks. Two-particle angular correlations originating from heavy-flavor particles allow such processes to be characterized.

The typical structure of a two-particle angular correlation distribution of high transverse-momentum (\(p_{\mathrm{{T}}}\)) trigger particles with associated charged particles features a “near-side” (NS) peak at \((\Delta \varphi , \Delta \eta ) = (0,0)\) and an “away-side” (AS) peak at \(\Delta \varphi = \pi \), extending over a wide pseudorapidity range. The NS peak is mainly induced by particles emerging from the fragmentation of the same parton that produced the trigger particle. The AS peak is related to the fragmentation of the other parton produced in the hard scattering. Here, \(\Delta \eta \) is the difference in pseudorapidity between the trigger and associated particles. The peaks lie on top of an approximately flat continuum extending over the full \((\Delta \varphi , \Delta \eta )\) range [47]. At leading order (LO) accuracy in QCD, heavy quark–antiquark pairs are produced back-to-back in azimuth [48]. At next-to-leading order (NLO), the correlation shapes can significantly differ from such a topology [48, 49]. Gluon radiation of heavy quarks can smear the back-to-back topology and broaden the near- and away-side peaks. In the gluon splitting process, the two heavy quarks can be produced with a small opening angle, depending on the \(p_{\mathrm{{T}}}\) of the gluon and the mass of the produced quark, generating two sprays of hadrons that can partially overlap, leading to a broader near-side peak. In the flavor excitation process [49], the heavy-quark pairs can be significantly separated in rapidity, and the hadrons from the opposite quark with respect to the trigger particle induce a nearly flat contribution to the \(\Delta \varphi \) distribution. The correlation measurements provide insight into heavy-flavor jet properties at low transverse momentum. By varying the \(p_{\mathrm{{T}}}\) interval of the trigger and associated particles, the correlation measurements allow the details of jet fragmentation to be studied, such as the jet angular profile and the momentum distribution of the particles produced in the fragmentation of the hard parton.

The azimuthal correlation distributions of prompt D mesons and charged particles were measured by the ALICE Collaboration in pp collisions at \(\sqrt{s} = 5.02, 7,\) and 13 TeV for \(p_{\textrm{T}}^{\textrm{D}}\) of the D mesons up to 36 \(\textrm{GeV}/c\)  and associated charged particles up to \(p_{\mathrm{{T}}}^{\mathrm{{assoc}}} =\) 3 \(\textrm{GeV}/c\) [47, 50, 51]. The measurements were compared with Monte Carlo (MC) simulations with different event generators, like PYTHIA [52,53,54], HERWIG [55, 56], EPOS [57, 58], and POWHEG coupled with PYTHIA8 for the parton shower and hadronization (POWHEG+PYTHIA8) [59, 60]. A substantial difference among the generators was observed, with PYTHIA8 and POWHEG+PYTHIA8 providing the best description of the measured observables. These differences can be ascribed to the specific implementation of features such as hard-parton scattering matrix elements, parton showering, hadronization algorithm, and underlying event generation, affecting the correlation functions of heavy-flavor hadrons and charged particles. Measuring the correlation distribution between heavy-flavor decay electrons and charged particles grants a substantially larger sample of correlation pairs, compared to measurements of D mesons and charged particle azimuthal correlations [47, 50]. This allows a significant extension of the \(p_{\textrm{T}}^{\textrm{assoc}}\) range of associated particles and can provide a more complete picture of the heavy quark fragmentation. In addition, electrons originating from beauty-hadron decays (b \(\rightarrow \) (c \(\rightarrow \)) e) dominate the heavy-flavor hadron decay electron spectrum (\(> 50\%\)) at high \(p_{\textrm{T}}^{\textrm{e}}\) (\(>5\) \(\textrm{GeV}/c\)) [61]. Hence, probing large enough trigger electron transverse momenta enables the study of the correlation function of particles originating from beauty-hadron decays, and provides information on the different correlation structures for charm and beauty quarks. This additional information can be used to further constrain the MC simulations. These advantages come at the price of an additional smearing introduced in the correlation function, due to the non-zero angle between the trigger electron direction and the direction of the parent heavy-flavor hadron before its decay. The momentum of the electron could also be further away from the quark momentum as compared to that of the parent hadron due to its decay kinematics.

In proton–nucleus (p–A) collisions, several cold nuclear-matter effects can influence the production, fragmentation, and hadronization of heavy quarks [5]. In the initial state, the parton distribution functions (PDFs) are modified in bound nucleons as compared to free nucleons. This feature is described by phenomenological parameterizations referred to as nuclear PDFs (nPDFs) [62,63,64]. When the production process is dominated by gluons at low Bjorken-x, the nucleus can be described by the Color-Glass Condensate (CGC) effective theory as a coherent and saturated gluonic system [65,66,67,68]. The CGC predicts momentum correlations in the initial state, that would impact the angular correlations of the produced heavy-quark pairs. Partons can also undergo multiple elastic, inelastic, and coherent scatterings, due to the presence of the nucleus in the initial state [69, 70] and to possible parton interactions in the high-density environment in the final state, particularly in collisions with large charged-particle multiplicity. These effects can be studied by measuring modifications in the angular shape or in the associated-particle peak yields of the angular correlation distributions of heavy-flavor particles with charged hadrons [47, 50]. Measurements of azimuthal correlations of prompt D mesons and charged hadrons in p–Pb collisions by the ALICE collaboration [47, 50], showed that the near- and away-side peaks of the correlation distribution are consistent with those measured in pp collisions in the same kinematic region. Employing heavy-flavor decay electrons as trigger particles in place of prompt D mesons allows studying the impact of cold-nuclear-matter effects for a wider associated particle \(p_{\textrm{T}}^{\textrm{assoc}}\) range, as well as to investigate their impact on the beauty-quark fragmentation and hadronization.

In heavy-ion collisions, a strongly-interacting matter consisting of deconfined quarks and gluons, the quark−gluon plasma (QGP), is produced [5, 71,72,73,74,75]. In the presence of the QGP, high-\(p_{\mathrm{{T}}}\) partons lose energy via medium-induced gluon radiation and collisions with the medium constituents [76,77,78,79,80,81]. These interactions cause a modification of the heavy-quark fragmentation and induce a broadening of the emerging jets and a softening of their constituents [82, 83]. Two-particle angular correlations have been extensively used to search for remnants of the radiated energy and to probe the medium response to the high-\(p_{\mathrm{{T}}}\) parton. The recent measurement of angular correlations between D mesons and charged particles in Au–Au collisions by the STAR Collaboration [84], shows a significant modification of the near-side peak width and associated yield, which increases from peripheral to central collisions. Measurements of angular correlations between electrons from heavy-flavor hadron decays and charged particles by the PHENIX Collaboration show modifications of the away-side peak yield and width in Au–Au collisions compared to pp collisions [85]. For future studies of heavy-flavor hadron correlations in heavy-ion collisions at the LHC, similar measurements in pp and p–Pb collisions are crucial to serve as reference [86].

In this article, ALICE measurements of the azimuthal correlations between electrons from heavy-flavor hadron decays with associated charged particles in pp collisions at center-of-mass energy \(\sqrt{s}=5.02\) TeV and p–Pb collisions at center-of-mass energy per nucleon–nucleon collision \(\sqrt{s_{\textrm{NN}}} = 5.02\) TeV are reported. The correlation distributions are measured for trigger electrons originating from heavy-flavor hadron decays in the \(p_{\textrm{T}}^{\textrm{e}}\) range \(4< p_{\textrm{T}}^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  and associated charged particles in the range \(1< p_{\textrm{T}}^{\textrm{assoc}} < 7\) \(\textrm{GeV}/c\), the latter granting a significantly higher \(p_{\textrm{T}}^{\textrm{assoc}}\) reach compared to previously published correlation measurements of D mesons with charged particles [47, 50]. The correlation distributions for trigger electron \(p_{\textrm{T}}^{\textrm{e}}\) in the range \(4< p_{\textrm{T}}^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< p_{\textrm{T}}^{\textrm{e}} < 16\) \(\textrm{GeV}/c\)  are also measured in order to study correlation shapes in kinematic ranges where the electrons are dominantly produced by charm- and beauty-hadron decays, respectively.

The article is organized as follows - in Sect. 2, the ALICE apparatus, its main detectors used in the analyses, and the data samples are reported. The complete analysis procedure is described in Sect. 3. The systematic uncertainties associated with the measurements are discussed in Sect. 4. The analysis results are presented and discussed in Sect. 5. The article is briefly summarized in Sect. 6.

2 Experimental apparatus and data samples

The ALICE apparatus consists of a central barrel, covering the pseudorapidity region \(|\eta | < 0.9\), a muon spectrometer with \(-4< \eta < -2.5\) coverage, and forward- and backward-pseudorapidity detectors employed for triggering, background rejection, and event characterization. A complete description of the detector and an overview of its performance are presented in Refs. [87, 88]. The central-barrel detectors used in the analysis are the Inner Tracking System (ITS), the Time Projection Chamber (TPC), and the electromagnetic calorimeters (EMCal and DCal). They are embedded in a large solenoidal magnet that provides a maximum magnetic field of \(B = 0.5\) T parallel to the beam direction. The ITS [89] consists of six layers of silicon detectors, with the innermost two composed of Silicon Pixel Detectors (SPD). The ITS was used to reconstruct the primary vertex and the charged particle tracks. The TPC [90] is a gaseous chamber capable of three-dimensional reconstruction of charged-particle tracks, and is the main tracking detector of the central barrel. Moreover, it enables charged-particle identification via the measurement of the particle specific energy loss (dE/dx) in the detector gas. The EMCal and DCal detectors [91, 92] are shashlik-type sampling calorimeters consisting of alternate layers of lead absorber and scintillator material. The EMCal covers ranges of \(|\eta | < 0.7\) in pseudorapidity and \(\Delta \varphi = 107^{\circ }\) (\( 80^{\circ }< \varphi < 187^{\circ }\)) in azimuth. The DCal is located azimuthally opposite the EMCal, with a coverage of \(0.22< |\eta | < 0.7\) and \(\Delta \varphi = 60^{\circ }\) (\( 260^{\circ }< \varphi < 320^{\circ }\)) and \(|\eta | < 0.7\) and \(\Delta \varphi = 7^{\circ }\) (\( 320^{\circ }< \varphi < 327^{\circ }\)). For the remaining part of this article, EMCal and DCal will be together referred to as EMCal, as they are part of the same detector system, used for electron identification. Two scintillator arrays, the V0 detector [93], placed on each side of the interaction point (with pseudorapidity coverage \(2.8< \eta < 5.1\) and \(-3.7< \eta < -1.7\)) were utilized for triggering and offline rejection of beam-induced background events. The minimum bias trigger was defined requiring coincident signals in both scintillator arrays of the V0 detector. In p–Pb collisions, the contamination from beam-induced background interactions and electromagnetic interactions was further removed with the information of the Zero Degree Calorimeters (ZDC) [94], located along the beam line at 112.5 m on both sides of the interaction point. A T0 detector [95], composed of two arrays of quartz Cherenkov counters, covering an acceptance of \(4.6< \eta < 4.9\) and \(-3.3< \eta < -3.0\), was employed to determine the luminosity together with the V0 detector.

The results presented in this paper were obtained using minimum bias triggered data recorded with the ALICE detectors during the LHC Run 2 from pp collisions at \(\sqrt{s}=5.02\) TeV and from p–Pb collisions at \(\sqrt{s_{\mathrm{{NN}}}}=5.02\) TeV. Pile-up events containing two or more primary vertices were rejected using an algorithm based on the detection of multiple vertices reconstructed from track segments in the SPD. In order to obtain a uniform acceptance of the detectors, only events with a reconstructed primary vertex within \(\pm 10\) cm from the center of the detector along the beam line were considered for both pp and p–Pb collisions. The number of selected pp and p–Pb events are about 800 M and 546 M, respectively, corresponding to integrated luminosities of \((16.63 \pm 0.32)\) nb\(^{-1}\) [96] and \((250 \pm 10)\) \(\mathrm {\mu b^{-1}}\) [97].

3 Analysis overview

The measurements of two-particle azimuthal correlations between electrons from heavy-flavor hadron decays (trigger) and charged (associated) particles were obtained from the correlation distributions of all identified electrons after subtracting the contributions which do not originate from heavy-flavor hadron decays. Effects from the limited two-particle acceptance and detector inhomogeneities were corrected using the event-mixing technique. The per-trigger correlation distributions were corrected for the associated-particle reconstruction efficiency. They were not corrected for the trigger-electron efficiency, as the efficiency was found to be \({ p}_{\textrm{T}}\) independent, and the correction factor would cancel with the per-trigger normalization. The properties of the correlation distribution in \(\Delta \varphi \), peak yields and widths, were obtained by applying a fit to the corrected \(\Delta \varphi \) distribution. A detailed description of the above mentioned analysis procedures is provided in the following sections. The analysis technique is the same in both pp and p–Pb measurements (unless specified otherwise in the text). Throughout this paper, the term “electron” refers to both electrons and positrons.

3.1 Electron identification and associated-particle reconstruction

Electrons with transverse momentum in the interval \(4< p_{\textrm{T}}^{\textrm{e}} < 16\) \(\textrm{GeV}/c\)  and \(|\eta | < \) 0.6 were selected using similar criteria as those discussed in Ref. [6]. Tracks were required to have at least one hit in any of the two SPD layers in order to reduce the contamination of electrons from photon conversions in the detector material. In order to reject secondary electrons [98], produced in interactions with the detector material or from weak decays of long-lived particles, the tracks were required to have a distance of closest approach to the primary vertex of less than 1 cm along the beam axis and 0.5 cm in the transverse plane. To ensure the selection of high-quality tracks, electron tracks were required to have a minimum of 70 crossed pad rows in the TPC (out of 159) and a minimum fraction of 0.8 of found space points relative to the maximum value, driven by the track direction [99]. The particle identification employed a selection on dE/dx inside the TPC and on the energy deposited in the EMCal detector. The discriminant variable used for the TPC detector is the deviation of dE/dx from the parameterized Bethe–Bloch expectation value for electrons [100], expressed in terms of dE/dx resolution, \(n\sigma ^{\mathrm{{TPC}}}_{\mathrm{{e}}}\). An asymmetric selection of \(-1< n\sigma ^{\mathrm{{TPC}}}_{\mathrm{{e}}} < 3\) was applied as the background contamination is higher for negative \(n\sigma ^{\mathrm{{TPC}}}_{\mathrm{{e}}}\). Additionally, electrons were identified and separated from hadrons using the E/p information from the EMCal detector, where E is the energy of the EMCal cluster (deposited by the particle while crossing the detector) [6, 101], and p is the momentum of the track measured by the TPC, along with a condition on the elliptical shape of the EMCal cluster, \(\sigma _{\mathrm{{long}}}^{2}\) [101]. The electron sample was obtained by selecting candidates with \(0.8< E/p < 1.2\), as expected for electrons, while hadrons have lower E/p values, and with \(0.02< \sigma _{\mathrm{{long}}}^{2} < 0.9\). The lower threshold on \(\sigma _{\mathrm{{long}}}^{2}\) removes contamination caused by neutrons hitting the readout electronics.

Associated particles were defined as all charged primary particles [98] with pseudorapidity \(|\eta | < 0.8\) and \({ p}_{\textrm{T}} > 1\) \(\textrm{GeV}/c\). Reconstructed tracks were required to have a minimum of 60 crossed pad rows in the TPC (out of 159) and a minimum fraction of 0.6 of found space points relative to the expected maximum considering the track position in the detector geometry [99]. Additional requirements on the distance of closest approach to the primary vertex of less than 1 cm along the beam axis and 0.5 cm in the transverse plane were applied. The associated particles were also required to have a \({ p}_{\textrm{T}}\) smaller than the trigger electron \({ p}_{\textrm{T}}\). This condition induces a kinematic bias for the regions where the trigger and associated \({ p}_{\textrm{T}}\) ranges overlap, that can be reproduced by simulations and model predictions.

3.2 Azimuthal correlation distribution and mixed-event correction

The two-dimensional correlation distribution as a function of azimuthal angle difference (\(\Delta \varphi = \varphi _{\mathrm{{e}}} - \varphi _{\mathrm{{ch}}}\)) and pseudorapidity difference (\(\Delta \eta = \eta _{\mathrm{{e}}} - \eta _{\mathrm{{ch}}}\)) between electron and charged particles, \(C(\Delta \varphi , \Delta \eta )\), was computed for the \({ p}_{\textrm{T}}\) interval \(4< p_{\textrm{T}}^{\textrm{e}} < 12\) \(\textrm{GeV}/c\), as well as the two intervals \(4< p_{\textrm{T}}^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< p_{\textrm{T}}^{\textrm{e}} < 16\) \(\textrm{GeV}/c\), and for five \({ p}_{\textrm{T}}\) intervals of associated particles between 1 and 7 \(\textrm{GeV}/c\)  (\(1< p_{\textrm{T}}^{\textrm{assoc}} < 2\) \(\textrm{GeV}/c\), \(2< p_{\textrm{T}}^{\textrm{assoc}} < 3\) \(\textrm{GeV}/c\), \(3< p_{\textrm{T}}^{\textrm{assoc}} < 4\) \(\textrm{GeV}/c\), \(4< p_{\textrm{T}}^{\textrm{assoc}} < 5\) \(\textrm{GeV}/c\), and \(5< p_{\textrm{T}}^{\textrm{assoc}} < 7\) \(\textrm{GeV}/c\)). For each kinematic interval, the correlation distributions were corrected for the limited pair acceptance and for the detector inhomogeneities using the event-mixing technique [102] as shown in Fig. 11 in Appendix A. The mixed-event correlation distribution, \(ME(\Delta \varphi , \Delta \eta )\), was obtained by correlating electrons in an event with charged particles from other events with similar multiplicity and primary-vertex position along the beam direction. The distribution obtained from the mixed events features a triangular-like shape as a function of \(\Delta \eta \), due to the limited \(\eta \) coverage of the detector, and is approximately flat as a function of \(\Delta \varphi \). Any non-flatness in \(\Delta \varphi \) would be due to \(\varphi \)-dependent detector inefficiencies and inhomogeneities. At \((\Delta \varphi , \Delta \eta ) \approx (0,0)\), the trigger and associated particle experience the same detector effects and the per-trigger correlation distribution is thus not affected. This property can be used to obtain the normalization factor, \(\beta \), for the mixed event distribution, defined as the average number of counts in the range \(-0.2< \Delta \varphi < 0.2\) and \(-0.07< \Delta \eta < 0.07\).

The mixed-event corrected correlation distribution, d\(^2N/(\)d\(\Delta \eta \)d\(\Delta \varphi \)), labeled as \(S(\Delta \eta , \Delta \varphi )\), was obtained as the ratio of the correlation distribution from the same event to the mixed event distribution, scaled by \(\beta \), i.e.,

$$\begin{aligned} \frac{\textrm{d}^{2}N}{\textrm{d}\Delta \eta \textrm{d}\Delta \varphi } \equiv S(\Delta \eta , \Delta \varphi ) = \beta \times \frac{C (\Delta \eta , \Delta \varphi )}{ME(\Delta \eta , \Delta \varphi )}. \end{aligned}$$
(1)

The two-dimensional correlation distribution was subject to significant statistical fluctuations, due to the limited size of the heavy-flavor decay electron sample, especially at large \(|\Delta \eta |\) values. To grant larger precision to the results, the mixed-event corrected azimuthal correlation distribution was integrated over in the range \(|\Delta \eta | < 1\) to obtain a one-dimensional \(S(\Delta \varphi )\) distribution.

3.3 Background subtraction

The hadron contamination in the selected electron sample was estimated by considering tracks identified as hadrons using \(n\sigma ^{\mathrm{{TPC}}}_{\mathrm{{e}}}<-3.5\). The E/p distribution of hadrons was scaled to match the electron-candidate E/p distribution in the interval \(0.3< E/p < 0.65\), away from the electron signal region, similar to the procedure discussed in Ref. [103]. The contamination from charged hadrons was estimated to be around \(1\%\) at \({ p}_{\textrm{T}}\) = 4 \(\textrm{GeV}/c\)  increasing to about \(12\%\) at 16 \(\textrm{GeV}/c\)  in both pp and p–Pb collisions. The hadron contamination in the azimuthal distribution of the inclusive electron sample was obtained using the correlation distributions of trigger particles with \(n\sigma ^{\mathrm{{TPC}}}_{\mathrm{{e}}} < -3.5 \), which was scaled to match the estimated hadron contamination. It was then subtracted from the inclusive electron (InclE) correlation distribution.

The selected electrons are composed of signal electrons originating from heavy-flavor hadron decays (HFe), and background electrons. The main background source is constituted by Dalitz decays of neutral mesons (\(\pi ^0\) and \(\eta \)) and photon conversions in the detector material, which produce electron–positron pairs with low invariant mass, peaked around zero. The background electrons were identified using an invariant-mass technique [104, 105], where each selected electron was combined with oppositely-charged partner electrons, obtaining unlike-sign (ULS) pairs and calculating their invariant mass (\( m_{\mathrm{e^{+}}e^{-}}\)). The partner electrons were selected by applying similar but looser track-quality and particle-identification criteria than those used for selecting the signal electrons, in order to increase the efficiency of finding the partner [105, 106]. Electron–positron pairs from the background have a small invariant mass, while random combinations including heavy-flavor decay electrons forming a pair with other electrons gives a wider invariant-mass distribution. This combinatorial contribution was estimated from the invariant-mass distribution of like-sign electron (LS) pairs. The \(S(\Delta \varphi )\) distributions of electrons composing ULS and LS pairs, \(S(\Delta \varphi )^{\mathrm{{ULS}}}\) and \(S(\Delta \varphi )^{\mathrm{{LS}}}\), respectively, were obtained. The background contribution was then evaluated by subtracting the LS distribution from the ULS distribution in the invariant mass region \(m_{\mathrm{e^{+}}e^{-}} < 0.14\) \(\textrm{GeV}/c\). The efficiency of finding the partner electron, referred to as the tagging efficiency (\(\varepsilon _{\textrm{tag}}\)) from here on, was estimated using MC simulations. In the pp and p–Pb analyses, the MC sample was obtained using PYTHIA 6.4.25 event generator [52], with the Perugia 2011 tune [107], and HIJING 1.36 [108] generators, respectively. They will be referred to as PYTHIA6 and HIJING in the following. The generated particles in all MC samples were propagated through the ALICE apparatus using GEANT 3.21.11 [109]. In order to increase the statistical precision of the tagging efficiency, \(\pi ^0\) and \(\eta \) meson samples with a flat \({ p}_{\textrm{T}}\) shape, generated with PYTHIA6, were embedded in the simulated events. The biased \({ p}_{\textrm{T}}\) shape was corrected by applying a weight to reproduce the measured \({ p}_{\textrm{T}}\)  spectra as described in [103, 110]. The tagging efficiency for pp (p–Pb) collisions was about \(74\%\) (\(75\%\)) at \({ p}_{\textrm{T}}\) = 4 \(\textrm{GeV}/c\), increasing to about \(79\%\) (\(77\%\)) for \({ p}_{\textrm{T}} > 7\) \(\textrm{GeV}/c\). The \(\Delta \varphi \) correlation distribution of background electrons was corrected by the tagging efficiency and subtracted from the inclusive electron distribution, that was already corrected for the hadron contamination, to obtain the azimuthal distribution of electrons from heavy-flavor hadron decays (\(S(\Delta \varphi )^{\mathrm{{HFe}}}\)),

$$\begin{aligned} S(\Delta \varphi )^{\mathrm{{HFe}}} = S(\Delta \varphi )^{\mathrm{{InclE}}} - \frac{1}{\varepsilon _{\mathrm{{tag}}}} [S(\Delta \varphi )^{\mathrm{{ULS}}} - S(\Delta \varphi )^{\mathrm{{LS}}}] .\nonumber \\ \end{aligned}$$
(2)

Contributions from other sources, such as decays of \(J/\psi \) and kaons, are negligible in the \({ p}_{\textrm{T}}\)  ranges considered in this analysis [104].

The azimuthal correlation distribution of electrons from heavy-flavor hadron decays and charged particles has to be corrected for the inefficiencies in the reconstruction of the associated particles and for the contamination of secondary particles in the associated particle sample. The reconstruction efficiency for charged primary particles was obtained using a different MC sample without any embedded particles using PYTHIA6 [52] and HIJING [108] generators for pp and p–Pb collisions, respectively. The efficiency obtained was in the range 86–90% (85–92%) in the \(1< { p}_{\textrm{T}} < 7\) \(\textrm{GeV}/c\) interval  for pp (p–Pb) collisions.

The amount of contamination from secondary particles [98] was also estimated using the same MC simulations, and shows values in the range 2–4% in pp collisions and 4–6% in p–Pb collisions, for the \({ p}_{\textrm{T}}\) interval considered. The fully-corrected azimuthal-correlation distribution was divided by the number of electrons originating from heavy-flavor hadron decays (\(N_{\mathrm{(c, b) \rightarrow e}}\)), to obtain a per-trigger normalization, where \(N_{\mathrm{(c, b) \rightarrow e}}\) is expressed as

$$\begin{aligned} N^{\mathrm{(c, b) \rightarrow e}} = N_{\mathrm{{InclE}}} - \frac{1}{\varepsilon _{\mathrm{{tag}}}} [N_{\mathrm{{ULS}}} - N_{\mathrm{{LS}}}] . \end{aligned}$$
(3)

3.4 Characterization of the azimuthal distribution

In order to quantify the properties of the measured azimuthal correlation, the following fit function was used

$$\begin{aligned} f(\Delta \varphi ) = b + Y_{\textrm{NS}}\frac{e^{\kappa _{\textrm{NS}}\cos {(\Delta \varphi )}}}{2{\pi }I_{0}(\kappa _{\textrm{NS}})} + Y_{\textrm{AS}}\frac{e^{\kappa _{\textrm{AS}}\cos {(\Delta \varphi - \pi )}}}{2{\pi }I_{0}(\kappa _{\textrm{AS}})}. \end{aligned}$$
(4)

It is composed of two von Mises functions, to model circular data, describing the NS and AS peaks, and a constant term, b, describing the baseline, which is a free parameter. The terms \(\kappa _{NS}\) and \(\kappa _{AS}\) in the von Mises function are the measure of concentration of NS and AS peak, respectively, where \(1/\kappa \) is analogous to the variance \(\sigma ^{2}\), and \(I_{0}\) is the zeroth-order modified Bessel function evaluated at \(\kappa \). The parameters \(Y_{\textrm{NS}}\) and \(Y_{\textrm{AS}}\) represent the integral of the near- and away-side peaks, respectively. By symmetry considerations, the means of the NS and AS peaks are fixed to \(\Delta \varphi = 0\) and \(\Delta \varphi = \pi \), respectively. The baseline b represents the physical minimum of the \(\Delta \varphi \) distribution. The width (\(\sigma \)) of the peaks is given by

$$\begin{aligned} \sigma = \sqrt{-2\log \frac{I_{1}(\kappa )}{{I_{0}(\kappa )}}, } \end{aligned}$$
(5)

where \(I_{1}\) is the first-order modified Bessel function evaluated at \(\kappa \). The per-trigger yields of the NS and AS peaks were obtained by integrating the bin counts in the ranges \(-3\sigma _{\mathrm{{NS}}}< \Delta \varphi < 3\sigma _{\mathrm{{NS}}}\) and \(-3\sigma _{\mathrm{{AS}}}< \Delta \varphi -\pi < 3\sigma _{\mathrm{{AS}}}\), respectively, after subtracting the baseline value b from the distribution.

4 Systematic uncertainties

The \(\Delta \varphi \) correlation distribution and the per-trigger NS and AS yields and widths are affected by systematic uncertainties, related to the procedures used for electron-track selection, identification and subtraction of the hadron contamination, background-electron subtraction, associated-particle efficiency correction, mixed-event correction, and fitting routine applied to the correlation distribution. The uncertainties from each of these sources were estimated separately, by varying the selection criteria or by using an alternative approach to the one described in the previous section. For each variation, its effect on the NS and AS peak yields and widths was obtained by reevaluating these observables after fitting and subtracting the baseline of the resulting correlation distribution. The uncertainties were computed separately for each trigger electron and associated particle \({ p}_{\textrm{T}} \) range. The systematic uncertainties on the correlation distribution from associated-particle efficiency correction and mixed-event correction are considered as correlated in \(\Delta \varphi \). The remaining sources are considered as uncorrelated in \(\Delta \varphi \). A summary of the systematic uncertainties of the correlation distribution, NS and AS yields and widths for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  are reported in Tables 1 and 2 for pp and p–Pb collisions, respectively. The \(\Delta \varphi \) correlated and uncorrelated uncertainties are separately reported for the \(\Delta \varphi \) distribution, and the total uncertainty from all sources is reported for the peak yields and widths.

Possible biases related to the specific track quality selection for electrons used in the analysis were studied by varying the selection criteria [6]. An uncertainty of 1–2% on the correlation distribution was obtained as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\) in both collision systems. For the NS and AS yields, an uncertainty in the range 1–2% was estimated. The uncertainty from track selection on the NS and AS widths was found to be negligible.

The uncertainty due to the electron identification using the TPC and EMCAL signals was estimated by varying the selection criteria for \(n\sigma ^{\mathrm{{TPC}}}_{\mathrm{{e}}}\), E/p, and \(\sigma _{\mathrm{{long}}}^{2}\). The chosen variations change the efficiency by a maximum of \(\sim \) 20%. A total uncertainty from these sources of 2–5% was obtained for the correlation distribution as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) in pp and p–Pb collisions, for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\). The resulting uncertainties ranged between 2% and 6% for the NS and AS yields, and between 2% and 7% for the NS and AS widths.

Table 1 Systematic uncertainties of the correlation distribution, the peak yields, and their widths for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\) in pp collisions. The individual sources of systematic uncertainties depend on the associated particle \({ p}_{\textrm{T}} \). The values presented as a range correspond to the lowest and highest \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) interval. For the correlation distribution, the systematic uncertainty from the baseline estimation is given as absolute value, and the total uncertainties from correlated and uncorrelated sources are reported separately
Table 2 Systematic uncertainties of the correlation distribution, the peak yields, and their widths for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\) in p–Pb collisions. The individual sources of systematic uncertainties depend on the associated particle \({ p}_{\textrm{T}} \). The values presented as a range correspond to the lowest and highest \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) interval. The systematic uncertainty of the correlation distribution from the baseline estimation is given as absolute values. For the correlation distribution, the systematic uncertainty from the baseline estimation is given as absolute value, and the total uncertainties from correlated and uncorrelated sources are reported separately

The contribution from background electrons was estimated using the invariant-mass method. The systematic uncertainty of the procedure, mainly affecting the average tagging efficiency, was obtained by varying the selection criteria of the partner electron tracks, including the minimum \({ p}_{\textrm{T}}\) and the invariant-mass window of the electron–positron pairs. The variation affects the tagging efficiency by \(\sim \) 5%. A resulting systematic uncertainty of 1–2% was obtained as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) on the correlation distribution, the peak yields, and their widths for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\) in pp and p–Pb collisions.

The uncertainty related to the specific selection of associated particles was estimated by varying the charged track selection criteria, including a requirement of a hit in one of the two SPD layers of the ITS, and varying the selection on the distance of closest approach, which affects the secondary particle contamination. This uncertainty is considered correlated in \(\Delta \varphi \). For \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\), uncertainties of 1–2% and 2–3% were obtained for the correlation distribution in pp and p–Pb collisions, respectively. For NS and AS yields, an uncertainty of 1–3% and 1–4% was estimated for pp and p–Pb collisions, respectively. Uncertainties of less than 3% and 4% were obtained for the NS and AS widths in pp and p–Pb collisions, respectively.

Effects induced by the limited detector acceptance and its local inhomogeneities were corrected using the mixed-event technique. The normalization factor, \(\beta \), was varied by taking the integrated yield over the full \(\Delta \varphi \) range for \(|\Delta \eta | < 0.01\). A correlated uncertainty in \(\Delta \varphi \) of 1% was obtained for the correlation distribution and the peak yields in pp and p–Pb collisions, respectively. No uncertainty was assigned for the NS and AS widths.

The \(\Delta \varphi \) distribution can be affected in case of a non-zero \(v_2\) of HFe and charged particles. As there are no previous measurements of HFe \(v_2\) in minimum bias pp and p–Pb collisions, a conservative estimate was obtained using the measurements in 0–20% central p–Pb collisions in Ref. [111]. The inclusion of \(v_2\) has an impact of less than 1% on the baseline and peak yields, and does not modify the NS and AS widths.

Several checks were performed to study the stability of the fit to the correlation distributions. Alternative functions, i.e., a Gaussian and a generalized Gaussian, were used to fit the NS and AS peaks instead of the von Mises function. Alternative fits were also performed fixing the baseline value to the average of the points in the transverse region, defined as \(\pi /3< |\Delta \varphi | < \pi /2\), to study its stability given statistical fluctuations. In place of the default bin counting procedure, the NS and AS yields were obtained as the integral of the fit functions in the range \(-3\sigma _{\mathrm{{NS}}}< \Delta \varphi < 3\sigma _{\mathrm{{NS}}}\) and \(-3\sigma _{\mathrm{{AS}}}< \Delta \varphi -\pi < 3\sigma _{\mathrm{{AS}}}\). The overall systematic uncertainty was calculated by taking the maximum variation of the results. The uncertainty from the baseline estimation on the correlation distribution is quoted as absolute numbers affecting all \(\Delta \varphi \) bins by the same value. The uncertainty of the NS and AS yields and width varies in the range 4–9% and 10–11% for pp and p–Pb collisions, respectively, for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\).

Similar procedures were followed to estimate the systematic uncertainties from the above mentioned sources on the correlation distribution, NS and AS yields and widths for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\). The uncertainty values were found to be similar to those obtained for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  in both collision systems.

5 Results

5.1 Comparison of the results in pp and p–Pb collisions

The azimuthal-correlation distributions for \(|\Delta \eta |<1\) with trigger electron in the interval \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\) and for different associated particle \({ p}_{\textrm{T}}\) ranges together with their fit functions are shown in Fig. 1 (for selected \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges) for pp (top panels) and p–Pb (bottom panels) collisions. The correlated systematic uncertainties, from the associated particle selection and mixed-event correction, are reported as text for each \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) interval. The baseline is shown by the horizontal green line. The absolute systematic uncertainty of the baseline estimation is shown as a solid box at \(\Delta \varphi \sim -2\) rad. The near- and away-side peaks are well described by the von Mises fit function in all \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges. While the baseline contribution is higher in p–Pb collisions (due to the larger charged-particle multiplicity), its absolute value reduces with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) in both pp and p–Pb collisions. As a large fraction of the baseline is from the underlying event processes, the pairs contributing to it are dominated by low \({ p}_{\textrm{T}}\) particles.

Fig. 1
figure 1

The azimuthal-correlation distribution for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  fitted with a constant function for the baseline (green line) and von Mises functions for AS and NS peaks (grey curves) for different associated \({ p}_{\textrm{T}}\) ranges in pp collisions at \(\sqrt{s}= 5.02\) TeV (top panels) and p–Pb collisions at \(\sqrt{s_{\textrm{NN}}}= 5.02\) TeV (bottom panels). The statistical (uncorrelated systematic) uncertainties are shown as vertical lines (empty boxes). The uncertainties of the baseline estimation are shown as solid boxes at \(\Delta \varphi \sim -2\) rad

To compare the NS and AS peaks of the \(\Delta \varphi \) correlation distribution between pp and p–Pb collisions, the baseline-subtracted distributions from the two collision systems are shown together in Fig. 2, for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  and for different \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges. It can be seen that the peak heights of the NS and AS decrease with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\). A tendency for a more pronounced collimation of the NS peak with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) is visible. The profile of the correlation peaks is consistent in pp and p–Pb collisions within the statistical and systematic uncertainties. This indicates that cold-nuclear matter effects do not impact heavy-quark fragmentation and hadronization in the measured \({ p}_{\textrm{T}}\) range, in minimum bias collisions. This observation is consistent with previous measurements of D-meson correlations with charged particles [47, 50].

Fig. 2
figure 2

Comparison of azimuthal-correlation distribution after baseline subtraction for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  and for different associated \({ p}_{\textrm{T}}\) ranges in pp collisions at \(\sqrt{s} = 5.02\) TeV and p–Pb collisions at \(\sqrt{s_{\mathrm{{NN}}}} = 5.02\) TeV. The statistical (uncorrelated systematic) uncertainties are shown as vertical lines (empty boxes). The uncertainties of the baseline estimation are shown as solid boxes at \(\Delta \varphi \sim -2\) rad

To perform a quantitative comparison of the correlation peaks between pp and p–Pb collisions, the per-trigger NS and AS peak yields (first row) and widths (third row) are shown in Fig. 3, superimposed for the two collision systems, as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\). The ratios between pp and p–Pb yields (second row) and widths (fourth row) are also shown in this figure. The systematic uncertainties on the ratio of the yields and widths were obtained by considering all sources except for the baseline estimation as uncorrelated between pp and p–Pb collisions. The partially correlated uncertainty of the baseline estimation, obtained by using different fit functions, was estimated on the ratio. The total uncertainty was obtained by taking the quadratic sum of the correlated and uncorrelated uncertainties. While the NS and AS yields decrease with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for both pp and p–Pb collisions, the measured yields are consistent within uncertainties between the two collision systems for all the \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges, as can be seen in the ratio panels of Fig. 3. The decrease in yields with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) can be understood considering that the heavy quarks have, on average, a hard fragmentation into heavy-flavor hadrons. As the remaining energy of heavy quarks is limited, it is far more likely that the associated particles accompanying the decay electron are preferentially produced at lower \({ p}_{\textrm{T}} \). The NS width values tend to decrease with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\), with a value of about 0.3 at \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}=\) 1 \(\textrm{GeV}/c\)  and narrowing to a value of roughly 0.15 at 6 \(\textrm{GeV}/c\), with a significance of about \(3\sigma \), for both pp and p–Pb collisions. The significance is calculated on the difference between the widths in the lowest and highest \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) intervals, taking into account both statistical and systematic uncertainties. The AS widths are independent of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\), and have a value of about 0.5. The NS peak distribution is closely connected to the fragmentation of the jet containing the trigger particle. The narrowing of the NS width with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) indicates that higher \({ p}_{\textrm{T}}\) particles tend to be closer to the jet-axis, whose direction can be approximated by the trigger electron. This is in turn related to higher \({ p}_{\textrm{T}} \) emissions from the heavy quark being more collinear to it. The AS peak exhibits a lower sensitivity to the fragmentation of a specific heavy quark, as it can contain particles produced via the fragmentation of heavy quarks originating from processes, including next-to-leading order, that are not azimuthally back-to-back. These processes may have different relative fractions for different \({ p}_{\textrm{T}} \) of heavy quarks. In the case of gluon splitting, the AS peak can also include particles originating from the recoil gluon, which are not directly associated with the heavy quarks produced in the event. Even in back-to-back processes, the correlation between the transverse momentum of the trigger electron and that of the opposite-side heavy quark, responsible for generating the AS peak through fragmentation, is significantly weaker than for the near-side peak. The NS and AS widths are similar in pp and p–Pb collisions, as can be seen in the ratio plots.

Fig. 3
figure 3

Comparison of near- and away-side per-trigger yields (first row) and widths (third row) as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  in pp collisions at \(\sqrt{s} = 5.02\) TeV and p–Pb collisions at \(\sqrt{s_{\mathrm{{NN}}}} = 5.02\) TeV. The ratios between pp and p–Pb yields and widths are shown in the second and fourth row, respectively. The statistical (systematic) uncertainties are shown as vertical lines (empty boxes)

5.2 Comparison with predictions from MC event generators

Fig. 4
figure 4

Comparison of the azimuthal-correlation distribution with model predictions after baseline subtraction for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  in different \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges in pp collisions at \(\sqrt{s} = 5.02\) TeV. The statistical (uncorrelated systematic) uncertainties are shown as vertical lines (empty boxes). The uncertainties of the baseline are shown as solid boxes near \(\Delta \varphi \sim 0\) rad

Fig. 5
figure 5

Comparison of the azimuthal-correlation distribution with model predictions after baseline subtraction for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  in different \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges in p–Pb collisions at \(\sqrt{s_{\mathrm{{NN}}}} = 5.02\) TeV. The statistical (uncorrelated systematic) uncertainties are shown as vertical lines (empty boxes). The uncertainties of the baseline are shown as solid boxes near \(\Delta \varphi \sim 0\) rad

Fig. 6
figure 6

Near- and away-side per-trigger yields (first row) and widths (third row) as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  compared with predictions from PYTHIA8 Monash tune and EPOS3 in pp collisions at \(\sqrt{s} = 5.02\) TeV. The ratios between model predictions and data are shown in the second and fourth row for the yields and widths, respectively. The statistical (systematic) uncertainties are shown as vertical lines (empty boxes)

Fig. 7
figure 7

Near- and away-side per-trigger yields (first row) and widths (third row) as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\)  compared with predictions from PYTHIA8 Angantyr and EPOS3 in p–Pb collisions at \(\sqrt{s_{\mathrm{{NN}}}} = 5.02\) TeV. The ratios between model predictions and data are shown in the second and fourth row for the yields and widths, respectively. The statistical (systematic) uncertainties are shown as vertical lines (empty boxes)

The near- and away-side peaks of the azimuthal-correlation distribution in pp and p–Pb collisions are compared with predictions from different MC event generators. This allows verifying the implementation of the processes of charm- and beauty-quark production, fragmentation, and hadronization, which have an impact on the observables studied in this paper. The models used to compare the measurement in pp collisions are PYTHIA8 with the Monash tune [44, 52, 53, 112] and EPOS 3.117 [57, 58]. The PYTHIA8 event generator is widely used in particle physics, as it provides an accurate description of high-energy collisions. It is capable of generating both hard and soft interactions, initial and final-state parton showers, particle fragmentation, and multi-partonic interactions. It also incorporates color reconnection mechanisms to rearrange color connections between quarks and gluons during hadronization. The prediction of these models for correlations of D mesons with charged particles can be found in Refs. [47, 50]. The p–Pb measurements are compared with PYTHIA8 Angantyr [113, 114] and EPOS 3.117 [57, 58] models. The Angantyr [113, 114] model is used to simulate ultra-relativistic p–Pb collisions with the PYTHIA8 event generator. As PYTHIA8 does not natively support collisions involving nuclei, this feature is implemented in the Angantyr model, which combines several nucleon–nucleon collisions to build a proton–nucleus (p–A) or nucleus–nucleus (A–A) collision. In this model, some modifications are made over the dynamics of pp collisions. The Angantyr model improves the inclusive definition of collision types of the FRITIOF model [115, 116]. In this model, a projectile nucleon can interact with several target nucleons where one primary collision looks like a typical pp non-diffractive (ND) collision. However, other target nucleons may also undergo ND collisions with the projectile. The Angantyr model treats secondary ND collisions as modified single-diffractive (SD) interactions. For every p–A or A–A collision, nucleons are distributed randomly inside a nucleus according to a Glauber formalism similar to the one described in Ref. [117]. This model is able to correctly reproduce final-state observables of heavy-ion collisions, i.e., multiplicity and \({ p}_{\textrm{T}}\) distributions [118]. As collectivity is not incorporated in this model, its predictions serve as a baseline for studying observables sensitive to collective behavior in p–A and A–A systems. For PYTHIA8 simulations, the correlation distributions for electrons from charm- and beauty-hadron decays are obtained separately, and summed after weighting their relative fractions based on FONLL calculations [30, 61, 119, 120].

The EPOS3 event generator is largely used for the description of ultra-relativistic heavy-ion collisions. It employs a core-corona description of the fireball produced in these collisions: in the “core”, its inner part, a quark–gluon plasma is formed, which follows a hydrodynamic behavior, while in the external regions of the “corona” the partons fragment and hadronize independently. A study of radial flow performed with the EPOS3 event generator in proton–proton collisions at \(\sqrt{s}\) = 7 TeV [121] has shown that the energy density reached in such collisions is large enough to grant the applicability of the hydrodynamic evolution to the core of the collision.

In the models, the azimuthal correlation function of trigger electrons from charm- and beauty-hadron decays with charged particles is evaluated using the same prescriptions applied for data analysis in terms of kinematic and particle-species selections. The peak properties of the correlation functions are obtained by following the same approach employed in data, i.e., by fitting the distributions with two von Mises functions and a constant term.

In Figs. 4 and 5, the baseline-subtracted azimuthal-correlation distribution measured in pp and p–Pb collisions, reflected in the \(0< \Delta \varphi < \pi \) range, is compared with predictions from PYTHIA8 and EPOS3 generators for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\) in three different \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges. The comparison for the remaining \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges is shown in Appendix B. From this qualitative comparison, both MC generators give a good overall description of the data in all the \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) intervals, even though the EPOS3 predictions show some deviation from the measured NS and AS peaks in the highest \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) interval. The peak yields and widths extracted from the measured distribution are also compared with model predictions in Figs. 6 and 7 for pp and p–Pb collisions, respectively. From here on, PYTHIA8/Angantyr will be used to refer to PYTHIA8 Monash simulations in pp collisions and PYTHIA8 Angantyr simulations in p–Pb collisions together. PYTHIA8/Angantyr simulations provide NS and AS yields decreasing with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) and are consistent with the data within statistical and systematic uncertainties. The NS widths simulated using PYTHIA8/Angantyr decrease with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\), which are consistent with the data in both collision systems. The AS widths show a slightly decreasing trend with \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) that is consistent with data within statistical and systematic uncertainties in both collision systems. The NS and AS yields predicted by the EPOS3 model qualitatively describe the data within statistical and systematic uncertainties in pp collisions. In p–Pb collisions, the NS yield is overestimated at high \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) while the AS yield is consistent with data within statistical and systematic uncertainties. The EPOS3 simulations overestimate the NS widths and underestimate the AS widths for all \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) ranges in pp and p–Pb collisions.

5.3 Dependence of the correlation distribution on the \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\)

Fig. 8
figure 8

Comparison of NS and AS per-trigger yields (first row) and widths (third row) for two \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) ranges \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\), as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) in pp collisions. The ratios between the \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\) and \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\) yields and widths are shown in the second and fourth rows, respectively. The data are compared with PYTHIA8 Monash and EPOS3 predictions. The statistical (systematic) uncertainties are shown as vertical lines (empty boxes)

Fig. 9
figure 9

Comparison of NS and AS per-trigger yields (first row) and widths (third row) for two \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) ranges \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\), as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) in p–Pb collisions. The ratios between the \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\) and \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\) yields and widths are shown in the second and fourth rows, respectively. The data are compared with PYTHIA8 Angantyr and EPOS3 predictions. The statistical (systematic) uncertainties are shown as vertical lines (empty boxes)

Fig. 10
figure 10

Comparison of PYTHIA8 Monash prediction for NS and AS per-trigger yields (first row) and widths (third row) in the two \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) ranges \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\) for electrons from charm- and beauty-hadron decays, as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) in pp collisions. The ratios to c, b \(\rightarrow \) e yields and widths are shown in the second and fourth rows, respectively. The statistical uncertainties are shown as vertical lines

The relative fractions of electrons produced by charm- and beauty-hadron decays have a strong \({ p}_{\textrm{T}}\) dependence [61]. The fraction of electrons from beauty-hadron decays at \({ p}_{\textrm{T}} ^{\textrm{e}}\) = 4 \(\textrm{GeV}/c\)  accounts for about 40% of the HFe yield, increasing to 60–70% for \({ p}_{\textrm{T}} ^{\textrm{e}} >8\) \(\textrm{GeV}/c\). A dependence of the correlation distribution on the flavor of the quark from which the trigger electron originates can be expected, due to the different fragmentation of charm and beauty quarks and different fraction of LO and NLO processes involved in their production. The correlation distributions for electrons from a given quark flavor can also have a trigger-particle \({ p}_{\textrm{T}}\)  dependence due to the different energy of the original parton, and different relative contribution of LO and NLO production processes for the hard scattering producing the parton. These effects are studied by measuring the correlation distributions for trigger electrons in the \({ p}_{\textrm{T}}\) ranges \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\), where the latter \({ p}_{\textrm{T}} ^{\textrm{e}}\) range is dominated by electrons from beauty-hadron decays. The azimuthal correlation distributions for these two \({ p}_{\textrm{T}} ^{\textrm{e}}\) ranges are presented in Appendix A. The NS and AS yields and widths for the two \({ p}_{\textrm{T}} ^{\textrm{e}}\) intervals are obtained following the same procedure described in Sect. 3.

The comparisons of the yields (first row) and widths (third row) for the two \({ p}_{\textrm{T}} ^{\textrm{e}}\) bins are shown in Figs. 8 and 9 for pp and p–Pb collisions, respectively. The per-trigger NS and AS yields are systematically higher for the \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\)  range compared to the values obtained for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\), for both pp and p–Pb collisions. The ratio between the \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\) and \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\) yields is shown in the second row of Figs. 8 and 9. It can be observed that the yield is higher for the higher \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) interval, and the ratio increases from 1.3 at low \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) to \(\sim 10\) in the highest \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) interval, for both pp and p–Pb collisions. This can be explained by considering that higher-\({ p}_{\textrm{T}}\) electrons are typically produced by more energetic heavy quarks, and the additional parton energy on average leads to a larger number of associated fragmentation particles.

While the NS width values decrease with \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\), they are similar for the two trigger electron \({ p}_{\textrm{T}} \) ranges. The AS widths are also observed to be similar for the two trigger electron \({ p}_{\textrm{T}} \) ranges and to have an almost flat trend with \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\). It should be noted that the kinematic bias induced due to the condition of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}} < { p}_{\textrm{T}} ^{\textrm{e}}\) affects the correlation distributions for the two trigger electron \({ p}_{\textrm{T}} \) ranges differently. While none of the correlation distributions for higher \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) interval are affected by the bias, the distributions for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\) and \(4< { p}_{\textrm{T}} ^{\mathrm{{assoc}}} < 7\) \(\textrm{GeV}/c\) would miss some associated particles because of the selection condition.

The NS and AS yields and widths of the correlation distributions as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for the two \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) ranges are compared with PYTHIA8/Angantyr and EPOS3 MC simulations for pp and p–Pb collisions. The PYTHIA8/Angantyr predictions describe the data within uncertainties for both \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) ranges. The NS and AS yields from EPOS3 are consistent with data for both \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) intervals. The trend of NS width from EPOS3 is slightly flatter as a function of \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) compared to that of the data. Similar to what was observed for \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\), the NS width is overestimated, while the AS width is underestimated compared to data for both \({ p}_{\textrm{T}} ^{\textrm{e}}\) ranges. The ratio of the yields and widths of the two \({ p}_{\textrm{T}} ^{\textrm{e}}\) ranges are well described by both MC event generators.

To understand the effect of the different charm and beauty fragmentation on the observed \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\)  dependence, the correlation distributions were obtained for electrons from charm- and beauty-hadron decays separately for the two \({ p}_{\textrm{T}} ^{\textrm{e}}\) intervals using PYTHIA8 MC simulations. The NS and AS yields and widths of the correlation distributions for electrons from charm- and beauty-hadron decays, and their ratios to the combined ones (HFe), are shown in Fig. 10. For both \({ p}_{\textrm{T}} ^{\textrm{e}}\) intervals, the NS yields for trigger electrons from beauty-hadron decays are lower than those from charm-hadron decays, by about 5% for the first \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) interval, with a tendency for an increased difference for larger \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\), about 40% for the last \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) range. This can be expected due to the harder fragmentation of beauty quarks to beauty hadrons compared to that of charm quarks, with less energy remaining for the production of other particles in the parton shower. This indicates that the yield increase at higher \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) observed in Figs. 8 and 9 is largely due to the higher energy of the initial heavy quark. The NS and AS widths of the correlation distributions decrease with increasing \({ p}_{\textrm{T}} ^{\textrm{e}}\) for both charm- and beauty-hadron decays, but the widths for electrons from beauty-hadron decays are wider than for electrons from charm-hadron decays for both \({ p}_{\textrm{T}} ^{\textrm{e}}\) intervals. These two opposing effects lead to similar widths for the two \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) intervals in Figs. 8 and 9.

6 Summary

Measurements of azimuthal-correlation functions of heavy-flavor hadron decay electrons with charged particles in pp and p–Pb collisions at \(\sqrt{s_{\mathrm{{NN}}}} = 5.02\) TeV have been reported. The correlation distributions were obtained for trigger electrons in the range \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 12\) \(\textrm{GeV}/c\), and for different associated particle \({ p}_{\textrm{T}}\) ranges between 1 and 7 \(\textrm{GeV}/c\). The azimuthal distributions were fitted with a constant and two von Mises functions in order to characterize the near- and away-side peaks.

The evolution of the near- and away-side peaks of the correlation functions in pp and p–Pb collisions is found to be similar in all the considered kinematic ranges. This suggests that the modification of the fragmentation and hadronization of heavy quarks due to cold-nuclear-matter effects is indistinguishable within the current precision of the measurements. The extracted near- and away-side per-trigger yields and widths in pp and p–Pb collisions are presented as a function of associated particle \({ p}_{\textrm{T}}\), which provide access to the momentum distributions of the particles produced in the fragmentation of the hard parton, and allow for a differential study of the jet angular profile. The per-trigger yields decrease with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) and are consistent between pp and p–Pb collisions. While the near-side width tends to decrease with increasing \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\), the away-side width does not show a pronounced trend with \({ p}_{\textrm{T}} ^{\mathrm{{assoc}}}\) for both collision systems. The \(\Delta \varphi \) distributions, per-trigger yields, and widths in pp and p–Pb collisions are compared with predictions from PYTHIA8 (with Monash tune for pp and using the Angantyr model for p–Pb collisions), and EPOS3 Monte Carlo event generators. The PYTHIA8 predictions provide the best description of the data for both yields and widths of the near- and away-side peaks. For the current implementation of the EPOS3 model, the yields are similar to those obtained from data, while the near- and away-side widths are overestimated and underestimated, respectively.

The relative fractions of electrons from charm- and beauty-hadron decays have a strong \({ p}_{\textrm{T}}\)  dependence. This feature was exploited by studying the correlation distribution for the kinematic regions, \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) \(\textrm{GeV}/c\)  and \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\), where the latter \({ p}_{\textrm{T}} ^{\mathrm{{e}}}\) range is dominated by beauty-hadron decays.

For both collision systems studied, the per-trigger yields are systematically larger for the \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\) range compared to the \(4< { p}_{\textrm{T}} ^{\textrm{e}} < 7\) interval due to the larger energy of the initial heavy quark, which allows for the production of more particles in the parton shower. This effect dominates over the increased beauty-origin contribution of the trigger electrons in the \(7< { p}_{\textrm{T}} ^{\textrm{e}} < 16\) \(\textrm{GeV}/c\) range, which according to PYTHIA8 studies are characterized by lower correlation peak yields than those of electrons originating from charm. The near- and away-side widths are observed to be similar for both trigger electron \({ p}_{\textrm{T}}\) ranges, for pp and p–Pb collisions. PYTHIA8 studies indicates that this is due to competing effects, where the larger boost of the initial heavy quark leads to a stronger collimation of the peaks with increasing \({ p}_{\textrm{T}} ^{\textrm{e}}\) for both charm- and beauty-origin contributions, compensating the broader peak widths for trigger electrons originating from beauty-hadron decays, whose contribution increases with \({ p}_{\textrm{T}} ^{\textrm{e}}\).

The reported results constitute a reference for future measurements in Pb–Pb collisions at the same center-of-mass energy. The study of the modifications of the correlation functions in Pb–Pb collisions in the presence of QGP can provide a deeper understanding of heavy-quark dynamics inside the hot QCD medium [86].